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Preface 


This book and disc is the second in a series of Computer Illustrated Texts 
(CITs) on topics in Numerical Analysis. The first, which was called A 
Simple Introduction to Numerical Analysis, covered topics such as the 
solution of algebraic equations, recurrence relations, linear equations, 
simple quadrature and the numerical solution of ordinary differential 
equations. In addition to the text each book is complemented by a 
disc of computer programs which are used as dynamic illustrations of 
the concepts discussed in the text and also remove much of the tedious 
calculation which is involved with repetitive processes. It is this integra- 
tion of a standard text with computer software which makes a CIT. This 
second volume extends the material of A Simple Introduction to Numer- 
ical Analysis to consider topics in approximation. These topics include 
interpolation, i.e. fitting a smooth curve through specified points in or- 
der to determine non-tabular values, the approximation of functions by 
polynomials and rational functions, least squares fitting of data, numer- 
ical quadrature and the solution of boundary value problems involving 
differential equations using finite element techniques. 

It is intended that the reader should have access to a suitable com- 
puter as he/she is reading the text and that at appropriate times the 
suggested programs are run. No detailed knowledge of computer pro- 
gramming is required though it is helpful if the reader is familar with 
the keyboard and how to input standard functions. However, full details 
of how to use the package are given in Chapter 1. 

This CIT is intended to be independent of A Simple Introduction 
to Numerical Analysis, however, at certain points in the text cross ref- 
erences are made to the original package. The reader will find it useful 
if this package is available but it is by no means essential. 


RDH 
DAQ 
May 1988 
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> Chapter 1 


> Getting Started 


> 1.1 Introduction 


This book is the second in a series of Computer Illustrated Texts on 
Numerical Methods; the first was called A Simple Introduction to Nu- 
merical Analysis. A Computer Illustrated Text (CIT) consists of a soft- 
ware package which is integrated with a textbook. It is intended that 
the reader should, ideally, use the software frequently whilst working 
through the book, as the textbook and software components have been 
designed to complement and reinforce each other. However, it is also 
intended that the textbook component can be used without access to 
a computer, and that the software can be used on its own, for exam- 
ple, to generate demonstration programs or to solve problems beyond 
the scope of this introductory text. It is the combination of text and 
programs which gives an ideal teaching environment. Each CIT is self- 
contained in that all the programs which are required are either included 
in the text or supplied on the disc. However, the nature of the subject 
means that no short introductory textbook can include all the material 
required and so explicit references are given to other sources that may be 
needed. No ordering in terms of difficulty is implied between this CIT 
and the first, but some of the exercises in the text can be performed 
more easily if the disc associated with A Simple Introduction to Numer- 
ical Analysis is available. This first chapter describes how to run the 
programs, input data, store results and display graphics. 

It is assumed that the reader has a basic knowledge of how to use a 
microcomputer but no knowledge of programming languages is needed, 
only an ability to enter data when required. 
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>1.2 The software: an example 


The software associated with this book was originally intended for the 
BBC microcomputer, but, other versions of the software will be avail- 
able, and all instructions in the main part of the book are machine 
independent. However, in this section the features of the software will 
be introduced using the BBC version of the program INTERP. You 
may prefer to skip now to Chapter 2 as it is always possible to refer 
back to this chapter. 

Users of machines other than the BBC microcomputer can usefully 
follow this section through, making only minor obvious variations. They 
will also find that their disc contains some introductory remarks which 
point out the most important of such variations. 


s Y ¢] * oy* 
Select FEL SPAS ETS 


eriz range: i 


Vert range: (-8.5 
H=5 Cursor at: 


Figure 1.2.1 Initial screen setting of INTERP. 


For BBC microcomputer users with the disc version of the soft- 
ware, loading of the program INTERP is accomplished by inserting 
the disc in drive 0, holding down SHIFT and then pressing and releas- 
ing BREAK. This should put up a menu on the screen giving a list of all 
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the programs which are available. (On a network system it will be nec- 
essary to ask for advice regarding the appropriate loading instructions.) 
On some screens you may find that either the top line or the bottom line 
is obscured; this can be overcome by selecting the menu option Screen 
adjustment and then following the screen instructions to obtain the best 
setting. Once the menu is visible use the cursor keys to select the re- 
quired program. Now select the program INTERP. The screen should 
clear and after a short delay it should look like figure 1.2.1. The screen 
has been split into three windows. The top window is used for dialogue, 
the bottom window contains current settings and the middle window is 
used either for graphic or tabular display. Initially, the option prompt 
should read 


Select data options 


where the data option is in inverse colours and a flashing cursor appears 
under the ‘d’ of data. The data option enables you to input, edit, 
plot and store data. Other options can be selected by pressing normal 
character keys: the SPACE bar is very convenient for this purpose. 
Alternatively, the key marked | can be used to cycle forwards and the 
key marked | to cycle backwards through the options available. An 
option is selected by pressing RETURN, but do not press RETURN 
just yet. 

Whenever the program is loaded a default set of data is supplied 
and plotted; the five default points are shown in the display window by 
crosses. These values can be tabulated and a method for doing so will 
be described later. In the bottom window the current ranges of plotting 
are given together with the number of data points (N). It is also possible 
to input data using the cursor, the current position of which is also given 
in the bottom window. It is a uniform convention, used throughout the 
software, that all data which may be changed by the user will have its 
existing setting left unchanged if the user presses RETURN or ENTER. 
Where options need to be selected, the alternative may be viewed by 
pressing the SPACE bar repeatedly until the desired option is showing, 
then selected by pressing RETURN or ENTER. 

The purpose of the program INTERP is to construct smooth curves 
through given data points. However, in this chapter we will only be in- 
terested in how data can be input, examined, edited, plotted and stored. 
Therefore, ensure that the Select option reads data options and then 
press RETURN to select this option. A new prompt should appear: 


Data: from keys 
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Cycle through the options available at this point by pressing the SPACE 
bar repeatedly. If at any time you wish to quit press ESCAPE which 
will take you back to the following prompt: 


Action: Continue 


Try it and see. Since we wish to continue, accept the default setting by 
pressing RETURN; the program should then advance to the Select data 
options prompt as above. Press RETURN again and we arrive back at 
the Data prompt. 

Select the option Data: from keys by pressing RETURN. Next we 
have the option of clearing any existing data points. For the moment 
let us just add to the existing points by using the default setting, i.e. by 
pressing RETURN. Additional points can either be input in tabular or 
graphical form. Pressing RETURN again will select the option 


Data input format: tabular 


and display the current points. A flashing cursor should now appear at 
the bottom of the table and instructions concerning the input of data are 
given in the top window. Try adding the point (0.95, 0.95) by typing 
in each coordinate followed by RETURN. If you make a mistake you 
can delete characters using the DELETE key before correcting them 
but each coordinate is terminated by RETURN. To quit press Q. This 
process will return the program to the option Data: from keys after 
sorting the tabular points into ascending order. Now press RETURN 
to accept the current option of supplying data from the keyboard and 
then press RETURN again to retain the existing points. The program 
now advances to the option Data input format: tabular. Press the 
SPACE bar to change this to Data input format: graphical and then 
press RETURN. All the current data points are displayed on the default 
graphical settings. These can be changed at this point using the Graph 
data option. Cycle through the options available, superimpose, clear - 
same axes and clear - new axes. For example, the last of these options 
allows the user to set up different sets of axes. When this option is 
shown press RETURN. The flashing cursor now appears in the bottom 
window and it is possible to alter the range of each axis in turn. For the 
moment accept the default settings by pressing RETURN each time the 
flashing cursor moves to a new position in the bottom window. When 
all the required values have been input the axes will be redrawn. 

We will next investigate the graphical editing options. Firstly, se- 
lect the option Add point. It is now possible to add additional points 
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either by (i) specifying the x and y coordinates or (ii) using the cursor. 
For (i) first press ‘X’ and a flashing cursor will appear in the bottom 
window as the x coordinate of the cursor. Try inputting a value 0.05 
followed by RETURN. Now press ‘Y’ and input a y coordinate of 0.0 
followed by RETURN. The flashing cursor in the middle window is now 
at the position (0.05, 0.0). In order to add this point to the table press 
RETURN again, this will be confirmed by an increase in N in the bot- 
tom window. To try out method (ii), first select the Add point option 
again by pressing RETURN. Additional points can now be added by 
moving the cursor with the arrow keys. The increment by which the 
cursor moves is controlled by the ‘<’ (increase) or ‘>’ (decrease) keys 
and is displayed in the top window. Note that the two methods can be 
combined; for example, you may set the x coordinate and then adjust 
it using the cursor control keys. The x and y values are open to change 
up to the moment when you press RETURN. 

The other graphical editing options are Change point, Delete point 
and Quit editor. The Change point option first asks you to select a point 
by using the cursor control keys; at this stage each keystroke makes the 
cursor jump to the next point in ascending order (right or up arrow), or 
descending order (left or down arrow). When the cursor flashes on the 
point you wish to change press RETURN, and the situation is exactly as 
it was with the Add point option. The Delete point option also begins 
by asking you to select a point; once you have selected it you are asked 
to confirm that you wish to delete the current point. The Quit editor 
option takes the program back to the option level from which the editor 
was invoked. 

At any stage of data input it is possible to view the current set of 
data points by escaping from the current option and selecting the option 


Data: tabulate 


An option to send any tabulation to a connected printer is available 
which reminds you to ensure that the printer is switched on. The option 


Data: plot &/or edit 


is available to display or modify any existing data set. 
It is possible to store data in a file using the option 


Data: store on file 
and retrieve suitable data using 


Data: from file 
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Wherever possible data files created by one program can be used as 
input for any other when it is appropriate, but if an inappropriate data 
file is requested then the program will display a suitable error message. 
Five simple data files are provided on the program disc!, as follows: 


DATA 
SIN 
SPLINE 
RAMP 
FE 


When loading a file you will need to type the whole file name at every 
occasion, however, you need not distinguish between upper and lower 
case. You will not be able to store any data on the program disc which 
is supplied because it is write protected and so you will need an unpro- 
tected disc of your own. You can give the data files any name which is 
consistent with the rules of the filing system you are using. Try loading 
the file DATA, examining it, editing it and then storing the result. 


> 1.3 The software: general remarks 


The previous section demonstrated how to supply simple numerical data 
and also the use of cyclical prompts. Although the actual prompts 
and options will vary from program to program the repeated use of 
the RETURN key will always accept the current default setting and 
consequently the first time a program is run results from the default 
parameters should agree with those given in the text. Throughout the 
programs the use of ESCAPE will always return the program to the 
Action prompt from which the program can be restarted or terminated 
as required. The ultimate panic button on the BBC model B is the 
BREAK key; this should only be pressed in emergencies. For MS-DOS 
versions BREAK and CTRL/C have the same effect as ESC. 

An alternative to numeric input is to supply data by means of ex- 
pressions. For example, at any point where the exponential number ‘e’ 
is required it is possible to type the string EXP(1). However, if the 
string supplied cannot be evaluated then a message such as ‘bad input’ 
or ‘error evaluating function’ will be given. Such an error will not nec- 
essarily be fatal to the program: a suitable error message will be given 
and usually the user will be invited to type in the item again. Some 

1With BBC discs data files have a prefix D.; MS-DOS files will have a .DAT 


extension. 
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errors cannot be detected until later in the calculation: these will also 
give rise to a suitable message, but it will act like an ESC in that the 
program will jump to the Action prompt. 


> 1.4 Screen dumps 


It is frequently useful to obtain a permanent copy of the screen dis- 
play for later examination. Programs running under MS-DOS may use 
the key labelled Print Screen, or an equivalent such as Prnt Sc, to 
obtain hard copies of the screen, provided that a suitable printer has 
been attached and has been selected by using the GRAPHICS command 
(consult your manual). For BBC microcomputer users, there is, unfortu- 
nately, no universal standard for printers or the software to drive them. 
Therefore, rather than supply a program which gives a direct screen-to- 
printer dump we have provided a facility for copying the screen image 
to memory, usually a disc file. The dump file can then be used to recre- 
ate the screen image and users can then provide their own software to 
make a hardcopy on any printer available. The details of this process 
will depend on the hardware available and the system being used. For 
the BBC microcomputer to which a disc drive is connected it is possible 
to obtain a screen dump as follows. Whenever a dump option is shown 
press RETURN at which point the top window should show 


Dump to filename: 


Make sure that the disc in the drive is not write-protected and then 
type in a suitable file name, which must satisfy the usual BBC filename 
conventions, followed by RETURN. 

If an Epson printer, or any other device which is software compatible, 
is to be used then the screen dump can be copied to the printer using the 
program PRNTDMP which is available from the main menu. Selecting 
this program gives a self-explanatory dialogue for the user to follow. 

One final advantage to this method of obtaining screen dumps is 
that it is possible to keep a permanent copy of screen images and recall 
them to the screen for future demonstration. For BBC users this can be 
carried out from within the program or else in intermediate command 
mode. In the latter case if the file DMPF were the name given to a 
dumpfile created in mode 4 then it can be recalled to the screen, but 
not printed, using the command 


MODE 4:*LOAD DMPF 
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> 1.5 Final remarks 


If you want to use numerical methods in more complicated programs, 
as opposed to learning about them, and you wish to write programs in 
BBC BASIC, then A Mathematical Toolkit: Numerical Routines with 
Applications in Engineering, Mathematics and the Sciences (Harding 
1986) contains subroutines which can be built into your own programs. 
For graphical output, Graphs and Charts on the BBC Microcomputer 
(Harding 1982) contains a library of subroutines for this purpose. It is 
intended that in the near future the CIT series will contain a Mathemat- 
tcal and Graphical Toolkit for MS-DOS machines. 


> Chapter 2 


> Interpolation 


‘Interpolation - the art of reading between the lines.’ 


> 2.1 Polynomial interpolation 


We shall begin our discussion of interpolation by considering the follow- 
ing problem. 


Example 2.1.1 

In an experiment the yield of a chemical species is recorded as a func- 
tion of the concentration of another reagent which is present. As it is 
impossible to vary the concentration of the reagent continuously, only a 
finite set of different values can be tabulated. Let us assume that the 
following results have been obtained. 


Table 2.1.1 Data for example 2.1.1. 


Concentration of reagent Yield 
0.10 0.09956 
0.30 0.39646 
0.50 0.58813 
0.70 0.77210 
0.90 0.89608 
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We can now ask questions like: 


1. Is it possible to estimate the yield when the reagent has a concen- 
tration of 0.51? 


2. How accurate is such an approximation going to be? 


3. Are the results which are produced sensitive to perturbations or 
errors in the tabular values? Clearly, we would like small changes 
in the data to produce small changes in the results. If this is the 
case then the problem is said to be well-conditioned. 


4. Is it possible to determine a reagent concentration which gives a 
prescribed yield, for example, find a concentration such that the 
yield is 0.5? (This is the inverse problem to (1)). 


The problem of attempting to estimate intermediate values from a 
given table is called interpolation and is the basis for the following sec- 
tions. 

In order to provide a background for understanding the methods 
described in this chapter we shall begin by considering the problem of 
obtaining non-tabular values in elementary ways. We assume that we 
are given a table of values (z;,y;), i = 0,...,n, and that we need to 
determine a non-tabular point (x,y). For example, let us suppose that 
we wish to find an approximate value which corresponds to z = 0.51 
in table 2.1.1. Perhaps the simplest way is to join together successive 
points as shown in figure 2.1.1. In particular, between x = 0.5 and 
xz = 0.7 we have the straight line shown in figure 2.1.2 and we can find 
the coordinates of the point C by linear interpolation. 


From figure 2.1.2 the point C has as its y coordinate y(0.5) + CD. 
However, by similar triangles, 


CD _ BE 
AD AE 
therefore, 
BE (0.77210 — 0.58813) 
CD = AD— = 0.01~—_———_—— = 0. 
AE (0.7 — 0.5) ee 


to five significant figures. (Only five significant figures are quoted here 
since the original tabular values are only given to this accuracy.) Hence, 
y(0.51) is approximately 0.58813 + 0.00920 = 0.59733. 
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Figure 2.1.1 A simple linear interpolant for table 2.1.1. 


Figure 2.1.2 Linear interpolation between two points. 
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Instead of using the similar triangle argument above we could have 
determined an explicit expression for the equation of the straight line 
joining the points A and B in figure 2.1.2. This line will be of the form 


y=co+ciz, 


and passes through the points (0.5, 0.58813) and (0.7,0.77210). There- 
fore, 
0.58813 = cp + 0.5¢1 


and 


0.77210 = co + 0.7}. 


Eliminating co gives cy = (0.77210 — 0.58813)/(0.7 — 0.5) = 0.91985, 
which is the gradient of the line. The intercept co is then found to be 
at 0.12821, hence the equation of the line joining the points A and B is 


y = 0.12821 + 0.919852, 


from which y(0.51) = 0.59733 as before. Clearly, this very naive ap- 
proach is simple to use but ignores most of the information given in the 
table. Furthermore, it fits an unnatural function to the data supplied 
since the fitted function has a strange derived function. Between each 
pair of consecutive data points the fitted function is linear and there- 
fore has constant slope. The derivative of the piecewise linear function 
shown in figure 2.1.1 will therefore be piecewise constant and look like 
the function shown in figure 2.1.3. In order to overcome this problem we 
could use higher order functions rather than linear ones. For example, 
through any three points we can fit a unique parabola. 


Exercise 2.1.1 (Analytical/Computational) Find the parabola through 
the points (0.3,0.39646), (0.5,0.58813) and (0.7,0.77210). (Hint: the 
general parabolic function is y = a+ bz +4 cz?, and the three coefficients 
can be determined by substituting the three points into this equation to 
produce three linear equations which must be solved to find values for 
a, b and c.) 


If the idea of fitting a polynomial through a number of points is unfa- 
miliar, then it is important to work through exercise 2.1.1. In particular, 
check that the resulting quadratic function takes the correct values at 
z=0.3,2=0.5 and z = 0.7. It is not difficult to extend this process to 
work for arbitrary polynomials, and it is to this problem that we next 
turn our attention. 
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Figure 2.1.3 The derived function for figure 2.1.1. 


A polynomial of degree n has n + 1 coefficients and so we can fit 
a unique polynomial of degree n through n + 1 distinct tabular points, 
Zi, 71 = 0,...,n. Let us assume that we are given a table of values 
(zi,yi), 7 = 0,...,n, and try to fit such a polynomial. Let P, be the 
interpolating polynomial 


P,,(z) = ep +12 + coe? +... +¢,2". 
Now we select the coefficients c;, 7 =0,...,n, so that the equations 
Peay = un pies PO ITs (24.1) 


are satisfied. This gives a system of n+ 1 linear equations in the n+ 1 
unknown coefficients cj, 7 = 0,...,n. These linear equations may then 
be solved using, for example, Gaussian elimination with partial pivoting. 
(See Harding and Quinney (1986).) 


Example 2.1.2 
For the five values in table 2.1.1 we can fit a unique polynomial of degree 
4, i.e. 


P4(z) =co+e12 + cox? + c3x? + cgr?. 


Evaluating this polynomial at the five data points in table 2.1.1 gives 
the following system of linear equations: 


14 


Tae Cees ge ae co 
ik Wek Gist WSy OR C} 
V 0,50) 0576 0 52a 05? co 
1 0:78 0:72) 05772 00 C3 
12 0:9' 0.92500, 028 0.9" C4 
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0.09956 
0.39646 

=| 0.58813 |. (2.1.2) 
0.77210 
0.89608 


Using Gaussian elimination, the solution of this system of equations is 


[—0.159796, 3.164656, —6.499406, 8.274375, ~3.901563]" 


which produces the interpolating quartic 


P(x) = —0.159796 + 3.1646562 — 6.4994062x? + 8.2743752° — 3.901563e7. 


(You can easily check these results using the program GAUSS from 
Harding and Quinney (1986).) Evaluating P4 at the given points pro- 
duces the values in table 2.1.2. (Later we shall provide an alternative 


simpler method.) 


Table 2.1.2 Polynomial interpolation applied to table 2.1.1. 


x y Pa(z) Pa(z) —y 
0.1 0.09956 0.0995598 2.4x 10-7 
0.3 0.39646 0.3964597 2.8x 107" 
0.5 0.58813 0.5881297 3.1x 10-7 
0.7 0.77210 0.7720996 3.9x 10-7 
0.9 0.89608 0.8960794 5.7x 10-7 


Notice that the values in the third column of table 2.1.2 are displayed 
to seven decimal places even though the data are only given to five. This 
is to emphasize the fact that whenever this process is carried out it will 
be prone to numerical errors. For this example, if we round the results 
in column 3 to five decimal places we recover the initial data. Whether 
this is always the case will be examined later. The function obtained, 
P4(x), is shown in figure 2.1.4. Notice that between the tabulated points 
the interpolating polynomial is quite smooth and so we are justified in 
evaluating it at non-tabular points. For example, we can determine a 
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Figure 2.1.4 Quartic polynomial for table 2.1.1. 


value for y(0.51) by P4(0.51) = 0.59733892. However, this is not always 
the case as the next example will demonstrate. Outside the range of the 
tabular points it would be unwise to use the interpolating polynomial to 
extrapolate the given values. This can be seen by examining the curve 
in figure 2.1.4 for values of x in excess of 0.9 where one might expect a 
monotonic increasing function. 


Example 2.1.3 

The interpolating polynomial which passes through the points (—2,e~*), 
(—1,e71), (0,1), (1,e7'), and (2,e7*) is shown in figure 2.1.5. These 
values are obtained from the function exp(—z”) which is always positive 
but as figure 2.1.5 shows the same is not true of the interpolating poly- 
nomial for values of z in the intervals (—2,—1) and (1,2)! Furthermore, 
as the number of points increases the interpolating polynomial behaves 
less like the function exp(—z”). This is called Runge’s problem and will 
be discussed in more detail in the next chapter. (See also exercises 2.1.4 
and 3.4.4.) 


Exercise 2.1.2 (Computational) The program INTERP is an interac- 
tive program which can be used to construct interpolating polynomials 
through a set of data points. The values in table 2.1.1 are set up as the 
default tabulated points. Therefore, load the program INTERP and 
when the screen shows 


Select data options 
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Figure 2.1.5 Runge’s problem: see example 2.1.3. 


press the space bar which will change the option to Select method. 
Press RETURN and the option prompt will change to 


Method: Construct interpolating polyn. 


You can cycle through the various methods by repeatedly pressing any 
character key except RETURN; now select the option above by pressing 
RETURN, which will calculate the required interpolating polynomial 
passing through the points in table 2.1.1. First of all, you are allowed 
the option of selecting the graphics parameters which are used to plot 
the resulting polynomial. For the moment select the option 


Graph data: superimpose 


The interpolating polynomial is now drawn in the graphics window and 
the actual polynomial is given in the bottom window. (The coefficients 
are limited to five significant figures by the space available to display 
them.) We can now evaluate the required interpolating polynomial using 
the option 


Option: Evaluate interpolant 


Entering the value for z = 0.51 produces the same value as in example 
2.1.2. All computation is carried out to full accuracy; only the coeffi- 
cients are displayed to five significant figures. 


Exercise 2.1.3 (Computational) Use the program INTERP and any 
four values from table 2.1.1 to determine a cubic interpolant and com- 
pare the estimate for y(0.51) with the value determined above. This 
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can be carried out by quitting the current option by pressing Q and 
then selecting the option Edit data. Press RETURN in order to ac- 
cept the current plotting parameters and then cycle through the edit 
options available by pressing any key but RETURN. (Recall that you 
can cycle backwards using the key marked |.) When the current op- 
tion is Delete point press RETURN and use the cursor keys to select 
the point you wish to remove and then press RETURN again. Cycle 
through the options now available until the option is Quit editor and 
then press RETURN. Once again cycle through the options available un- 
til the current option becomes Select method and press RETURN yet 
again. Accepting the option Method: Construct interpolating polyn., 
the interpolating polynomial through the remaining points will be deter- 
mined and plotted. The effect of deleting another point can be examined 
as follows. Cycle until the option Edit data is current and then press 
RETURN twice. This will take you into the data editor; the current 
edit option should be Add point. If this is the case press RETURN, if 
not cycle until it is and then press RETURN. The cursor in the graph- 
ics window should now be flashing at the point which was previously 
deleted; this point can be added back to the table simply by pressing 
RETURN. Other points can then be deleted to consider the effect of 
fewer points in the table. How does the selection of the four points 
influence the interpolated value? 


Exercise 2.1.4 (Computational) Use the program INTERP to obtain 
five equally spaced points in the interval (—5,5) from y = exp(—2”) 
and then determine a suitable interpolating polynomial (see example 
2.1.3). You will need to take the default prompt Select data option 
and then Data: from keys remembering to clear the existing data be- 
fore using the Data input format: tabular option. For example, the 
point (—5, exp(—25)) can be input as the tabular values X(1) = —5 and 
Y(1) =EXP(—25). When you have constructed the interpolating poly- 
nomial through the points you have selected, try adding other points and 
see how the results compare. What happens as the number of points 
increases /decreases? Why is it better to have an even number of points? 
Investigate what happens with non-equally spaced points. 

In general, given n +1 tabular points, (z;,yi), 7 = 0,...,n, we 
can proceed as above and determine the coefficients of the interpolat- 
ing polynomial P,(x) by solving a system of linear equations given by 
equation (2.1.1), or equivalently 


Ac=y, (2.1.3) 
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where A is an n+1 by n+1 matrix and c and y are vectors of di- 
mension n + 1. The columns and rows of A are numbered 0 to n and 
the components of A are aj; = 27,i1=0,...,n, j =0,...,n. Similarly 
yi = y(z:), 1 = 0,...,n, and ¢; are the coefficients of the interpolat- 
ing polynomial which is required. (See equations (2.1.2) for a simple 
example with n = 4.) This system of equations will have a unique solu- 
tion provided the matrix A is non-singular, therefore, we examine the 
possibility of it being singular or almost singular. In order to do this 
we consider the determinant of A; this matrix is called a Vandermonde 
matriz and it is possible to show that it has determinant 


Lie aG@es,): 
Cet cca tn ieee 


Exercise 2.1.5 (Analytical) Show that the determinant of the Vander- 
monde matrix of order 3, given by 


lean 
ere cl 
1 ope aa 


is (x1 = Lo)(x2 — Lo)(X2 = Zils 


The importance of the above expression for the determinant of the Van- 
dermonde matrix is that if the tabulation points are close together then 
the determinant may be small, i.e. the equations are nearly singular, and 
hence the solution may be inaccurate. This type of inaccuracy is usually 
called inherent ill-conditioning and is a fundamental part of the problem, 
i.e. independent of the method of solution. In order to demonstrate the 
problems which may arise when determining interpolating polynomials 
in such circumstances let us consider the following example. 


Example 2.1.4 (Computational) 

Let us begin by finding the quartic polynomial which passes through 
the points in table 2.1.3. Using the program INTERP produces the 
interpolating polynomial shown in figure 2.1.6. In order to consider how 
sensitive this problem is let us perturb the y value of the third entry 
in table 2.1.3 from 1.05 to 0.95 and see how this affects the computed 
solution. The resulting quartic polynomial is shown in figure 2.1.7. Com- 
pared with the original interpolation quartic it is clear that this problem 
is inherently ill-conditioned. 
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Table 2.1.3 


eat spent 5 
0.00 2.00 
0.95 1.00 
1.00 1.05 
1.05 1.00 
2.00 0.00 


Figure 2.1.6 The quartic interpolant for table 2.1.3. 


Figure 2.1.7 The effect of perturbing entry 3 in table 2.1.3. 
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Example 2.1.4 gave an illustration of a problem which was inherently 
ill-conditioned, i.e. a small change in the tabular values leads to a much 
larger change in the interpolating polynomial. This was caused by sev- 
eral of the tabular values being close together, resulting in a system of 
linear equations which were nearly singular. We can also investigate the 
possibility of errors being introduced in the interpolating polynomial by 
the method we have selected to solve this problem. Such perturbations 
are called induced errors. Ideally we would like small changes in the 
tabular values to produce small changes in the interpolating polynomial 
and if this is the case then the problem is said to be well conditioned. 
In order to demonstrate this idea consider the following example. 


Example 2.1.5 (Computational) 
Applying the program INTERP to construct an interpolating polyno- 
mial through the points in table 2.1.1 produces the polynomial 


P,(x) =—0.159796 + 3.1646562 — 6.499406x? + 8.2743752° — 3.901563z%. 


Now let us round the entries in table 2.1.1 to four decimal places and 
then construct a polynomial through the resulting points. This gives the 
polynomial 


—0.1597 + 3.16522 — 6.50442? + 8.283322 — 3.906327. 


The maximum relative change in any coefficient is less than 0.12% and 
so this problem is relatively well conditioned. 


Exercise 2.1.6 (Computational) Use the program INTERP to construct 
the interpolating polynomial through the tabular values in the file SIN. 
Investigate the sensitivity of the data provided by examining the effect 
of perturbing one or more of the values given. Is this problem inherently 
well conditioned? Is it susceptible to induced ill-conditioning? 


> 2.2 Lagrangian interpolation 


In order to construct a polynomial of degree n through n+ 1 distinct 
points using equation (2.1.1) it is necessary to solve a system of linear 
equations. We will now derive an alternative way of producing such a 
polynomial which does not require the solution of a system of possibly 
ill-conditioned equations. As there is only one polynomial of degree n 
which passes through n +1 points we would obtain the same polynomial 


LAGRANGIAN INTERPOLATION 7211 


as that obtained from equation (2.1.1) if exact arithmetic were used. 
(See exercise 2.2.4.) 


Let us begin by assuming that we have been given a set of data 
points, (zj,y;), 7 = 0,...,n, and then examine the function 
(z — 20)(z — 21)... (@ — 23-1) (2 — 2441)...(2 — Zn). 


Notice that there is no factor (x—2;). This function can also be written 


as the product 

[]( - 2)). 

jFi 
and is a polynomial of degree n which vanishes at each data point with 
the exception z = 2;. We now define the Lagrangian polynomial func- 
HONS; (a )yt = Uyes 2, as 
NPAC: — 23) 
TT] zi (2 — 2;)’ 
We see that each of these polynomials has degree n and that L;(x,) = 0 
if k #17 but is equal to 1 when k = 7. Now consider the function 


Li(x) = i=0,...,n. (2.2.1) 


Pian) =  wLi(e). (2.2.2) 


Since L;(x) is a polynomial of degree n for each 7, so is P(x) and by 
virtue of the way in which each of the functions L; are defined P(z) 
will take the values y; when x = 2;, 1 = 0,...,n. Therefore, P(x) and 
the interpolating polynomial P,(z) are identical since they are both 
polynomials of degree n which agree at n+ 1 points. 


Example 2.2.1 (Analytical) 
Let us consider the data points shown in table 2.1.1. The Lagrangian 
polynomial Lo is given as follows: 
(x — 0.3)(2 — 0.5)(x — 0.7) (x — 0.9) 
(0.1 — 0.3)(0.1 — 0.5)(0.1 — 0.7)(0.1 — 0.9) 
= 26.04167x* — 62.52° + 53.645832* — 19.3752 + 2.460937. 


We stress that in practice an explicit expression of such a polynomial 
would not be determined; it is easier to evaluate it in the form of equa- 
tion (2.2.2). You should check that Lo(zo) = 1, i.e. Lo(0.1) = 1, and 
that Lo(2;) = 0, i = 1,2,3,4, ie. L0(0.3) = Lo(0.5) = Lo(0.7) = 
Lo(0.9) = 0. However, because of the form of equation (2.2.2) it is beter 
to consider y;L;(x), for example, yoLo(x) is plotted in figure 2.2.1. 


Lo(x) = 
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Figure 2.2.1 The Lagrangian polynomial yoLo(z). 


Similarly, we can construct the polynomials L(x), Lo(x), L3(x) and 
L4(x). The required interpolating quartic polynomial is then given by 


Py(x) = 0.09956L9(zx) + 0.39646L; (x) + 0.58813L2(z) 
+0.77210L3(x) + 0.89608L4(z). 


To evaluate P,(0.51) we firstly find Zo(0.51) = 0.004052 from equa- 
tion (2.2.1) and then, similarly, we have that D,(0.51) = —0.03165, 
L2(0.51) = 0.99688, L3(0.51) = 0.03498 and L4(0.51) = —0.00426, 
hence P4(0.51) = 0.59734 as before. 


Exercise 2.2.1 (Computational) The program INTERP may be used to 
construct a suitable interpolating polynomial through the points given 
in table 2.1.1. Load INTERP, cycle through the options available with 
the Select prompt until it reads method and then press RETURN. 
Cycle again until the Method prompt is Lagrangian polynomials and 
press RETURN again. Now use the default Graph option which is su- 
perimpose and a flashing cursor will appear in the graphics window. 
Different points can be selected using the cursor keys and then pressing 
RETURN once again. For example, the Lagrangian polynomial Lo(z) 
can be drawn by simply pressing RETURN when the flashing cursor is 
positioned over the first tabular point. However, in order to make use 
of the form of equation (2.2.2) it is more convenient to use the functions 
yoLo(xz) and it is this function which is plotted in the graphics window 
and displayed in the bottom window. In a similar way the other four 
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Lagrangian polynomials can be drawn. 


Erercise 2.2.2 (Computational) Use the program INTERP to construct 
the Lagrangian polynomials through the data in file DATA. Sum these 
polynomials as in equation (2.2.2) and show that the result is equiva- 
lent to the polynomial given by direct construction of the interpolating 
polynomial. 


Erercise 2.2.3 (Computational) Use the program INTERP to determine 
suitable third- and fourth-order polynomial interpolants for the data 
given in table 2.1.1. (Hint: use the editor to delete points one at a 
time.) 


Exercise 2.2.4 (Analytical) Show that there is a unique polynomial of 
degree n which passes through the (n + 1) points (2;,y;),7 = 0,...,n 
(Hint: consider two polynomials p,(z) = ag + a,x +...+a,2" and 
Qn(z) = bo +bin+...+ 6,2", with a; # 6; for at least one i, each 
of which interpolates the points given. Now consider the non-trivial 
polynomial pn(z) — gn(x) at the points z; and hence show that a; = ;, 
toma Meee. 1) 


Exercise 2.2.5 (Analytical) Hermite interpolation. Lagrangian interpo- 
lation is based upon tabular values (z;,y;), 7 = 0,...,n; if in addition 
the gradients y; are supplied then higher order polynomials can be fitted 
to the data. Show that it is possible to fit a polynomial of degree 2n + 1 


in the form 
Peale) = yal x)yi + ya) 


where 


[1 — 2L;(2i)(x — xi)]Li(x)’, 
(x = x;)L;(z)?. 


a;(zr) 


Pi (zx) 


The result is called a Hermite interpolation polynomial. 


> 2.3. Divided differences 


As another alternative to the construction of an interpolating polynomial 
through the points (z;,yi), 1 = 0,...,n, using Lagrangian polynomials 
we can determine the interpolating polynomial of degree n as follows. 
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If we define the function po(z) = yo then, trivially, po(zo) = yo. Now 
consider the linear polynomial 


pi(x) = yo + ai(z — Z0), 


which satisfies pi(ro) = yo and then select a; = (y — yo)/(x1 — Zo) so 
that 


pi(z1) = yo + 41(21 — Zo) = 1- 
Next define 
p2(z) = ao + a1(x — 20) + a2(x — z0)(x — 21), 


then since aj = yo and ay = (y; — yo)/(21—Zo) we have that po(xo0) = yo 
and po(z,) = y;. Furthermore, 
po(x) = pi(x) + ao(z — 29)(x — 24). 
We still have az at our disposal, but if 
Ss (p2(z2) — pi(%2)) 
(x2 — £0)(x2 — 21) 


then po(x2) = yz, i.e. p(x) passes through (zo0, yo), (21, yi) and (22, y2). 
We can repeat this process by defining 


pj (2) = pja1 (2) + a;(@ 20) (@ 23) (2 — £o)...(2 — 2524), J = Vent 


and then selecting a; so that p;(z;) = y;. The resulting polynomial, 
Pn(z), passes through (2;,y;), 7 = 0,...,n, and is identical with that 
obtained by using the Lagrangian polynomials. 


Example 2.3.1 (Analytical) 
The interpolating polynomial through the first three values in table 2.1.1 
may be constructed as follows. Using the above method gives 

po(z) = y= 0.09956, 

Pi(z) = po(z) + a;(z — zo) = 0.09956 + a, (zx — 0.1). 


Select a, so that p;(xz,) = 0.39646, i.e. ay = Sa CELE which gives 


Pi(z) = 0.09956 + 1.4845(2 — 0.1). 
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Next we define p2(z) by 


po(x) = pi(x) + ao(x — 0.1)(z — 0.3) 
and select aj so that po(x2) = yo, i.e. po(0.5) = 0.58813, 
(p2(t2) — pi(22)) 
(x2 — 20)(t2 — 21) 
(0.58813 — 0.69336) _ 
(0.5 —0.1)(0.5-—0.3) — 


This gives the polynomial 
po(x) = 0.09956 + 1.4845(2 — 0.1) — 1.315375(2 — 0.1)(z — 0.3) 


ag = 


—1.315375. 


which passes through the first three entries in table 2.1.1. 

The process outline above can be automated in the following way. 
Given the table of values (z;, y;), i = 0,...,n, we define the first divided 
differences associated with these values by 
(yita — Yi) 

(vi+1 — Zi) 
Next we compute second divided differences as 
(f[vi+1, 2142] — f[zi, 2i41]) 
(zi42 — 2:) 
and so on. When the (k — 1)th divided differences have been obtained, 
1.e. 


here Li41] = 


Tony Fearn Cone 


fle eaeer eee) Cid EoMaAN SL tee ts Cit 2c CALE, 
we define the kth divided difference as 
(flariy- ty toy r+ Bleed [235 Fi +k —1)) 
Seceei enna Cs = eee S"Wm_—_—0—omms. 
f| +1 +k] (xs4% ates z;) 


These differences can be shown schematically as in table 2.3.1. 


We can then show that the polynomial p,(2x), given by 
Pn(z) = yotf[zo,21])(2 — 20) + f[z0, 21, 22](e — xo)(z — 21) 
ey fortieeetn lf 220) f 7 21)-2(2 — fn=1)3 


passes through the points (2;, y;), 7 = 0,...,n, and is therefore identical 
to the interpolating polynomial P,(z). (See exercise 2.2.4.) 
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Table 2.3.1 Table of divided differences. 


oY First divided Second divided Third divided 
difference difference difference 

Zo Yo 

f[zo,21]= Gia 
T1 Yi f[2o,20eo) eee 

f[z1,r2]= fa) f[xzo,x1,£2,23] 
£2 Y2 f(z1,22,05] = Heazslaferaal 

f[x2,r3]= (arva) 
T3 Y3 


Exercise 2.3.1 (Analytical) Show that the above expression is the inter- 
polating polynomial through the points (z;,y;),7=0,...,n. 


Exercise 2.3.2 (Analytical) Show that f[r1,22,...,2,] is invariant with 
respect to all permutations of the tabular points 21, 22,..., 2x. 


Table 2.3.2 Tableau of divided differences for table 2.1.1. 


First Second Third Fourth 
i y difference difference difference difference 
0.1 0.09956 
1.48450 
0.3 0.39646 —1.31538 
0.95835 2.031875 
0.5 0.58813 —0.09625 —3.9015625 
0.91985 —1.089375 
0.7 0.77210 —0.74988 
0.61990 
0.9 0.89608 


Example 2.3.2 (Calculator) 


To demonstrate the use of divided differences we construct the inter- 
polating polynomial through the five points in table 2.1.1. Firstly, we 
obtain the table of divided differences given by table 2.3.2, as illustrated 
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in table 2.3.1. The required quartic polynomial is then given by 


P(x) = 0.09956 + 1.48450(x — 0.1) — 1.31538(x — 0.1)(x — 0.3) 
+2.031875(x — 0.1)(x — 0.3)(2 — 0.5) 
—3.9015625(x — 0.1)(x — 0.3)(x — 0.5)(x — 0.7), 


where the coefficients are given by the first entry in each column of 
the tableau of divided differences. Expanding this expression gives the 
polynomial 


P4(x) = —0.159796 + 3.1646562 — 6.4994062x? + 8.2743752° — 3.9015632%, 


which agrees with the interpolating polynomial through these points as 
computed previously. 


Exercise 2.3.3 (Analytical) Use the method given in section 2.1 to deter- 
mine an interpolating polynomial of degree 2 which passes through the 
points (0,3), (1,6) and (2,11). Obtain the corresponding Lagrangian 
polynomials Lo, L; and Lz and hence show that the polynomial 


3Lo(x) + 6L1(x) + 11L2(z) 


is identical with the interpolating polynomial which was derived using 
equation (2.1.1). Construct a table of divided differences for these data 
points and hence construct a quadratic function which passes through 
the same points. Confirm that all three methods give the same polyno- 
mial. 


Exercise 2.3.4 (Programming) Write a simple program to read in a set 
of tabular points and construct a table of divided differences. Test your 
program on the data given in exercise 2.3.3. Extend the program to 
produce polynomial interpolants as described in example 2.3.1. 


> 2.4 Neville’s algorithm 


Although the previous sections demonstrated how it was possible to 
construct. a suitable interpolating polynomial explicitly it is clear that 
a considerable amount of effort is required in order to interpolate from 
the values given in a table. Furthermore, neither the method outlined in 
section 2.1 nor the Lagrangian method gives any indication of the order 
of polynomial required. What we seek is some means of using a subset 
of the tabulated points, m say, to construct a suitable interpolant and 
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then by adding an additional point construct a higher order polynomial. 
In this way if we can show that the interpolants produced using subsets 
of m points and m+1 points agree then there is no point in investigating 
higher order polynomials. Neville’s algorithm does precisely this. 

Let us suppose we select a subset of m points from a table of val- 
ues, (20, Yo), (£1, Yi), ---» (Lm—1) Ym-1)- Next we construct the linear 
functions 


(z — Zo)yi — (2 — 21) yo 


Poi(z) 


(1 — 20) , 
a} (2 — 21) yo — (z — F2)y1 
Pio(z) = beh ee ne =) ’ 


(x “Fr EA) Vind aay (x a Zn) Ure 
(ore a Tim=2) 


Dine yeaa e 


Notice that p; ;41(2) passes through the points (2; ,y;) and (2;41,yj+41), 
j =0,...,m—1. (These terms are sometimes called linear cross means.) 
Next we construct the linear cross means between successive linear terms, 
e.g. between poi(x) and pi2(x), and then evaluate the quadratic function 


(z — x0)pi2(x) — (x - %2)por(x) 
(r2 — £0) 


This function passes through (Zo, yo), (21, yi) and (z2,y2) and is de- 
noted by poo(xz). In this way we can recursively define a sequence of 
polynomials of increasing order by 


(x — 2; )pj41,4(2)—(2 — te); k-1(2) 


Pjk (x)= i ej Roe eee ite 
oe ee 
(2.4.1) 
This polynomial is of degree (k—j) and passes through the points (z;, y;) 
fori = j,...,k, and is therefore equivalent to the Lagrangian polynomial 


of the same degree (see exercise 2.2.4). 
Exercise 2.4.1 (Analytical) Show that the quadratic polynomial po3(z) 
satisfies the following: po3(%o) = yo, pos(@1) = Yi, Pos(Z2) = Yy- 

In practice if we are given a set of data and require an estimate at 


a non-tabular point 2* then, rather than evaluating the general polyno- 
mial, we only evaluate (2.4.1) at the point z*. Furthermore, it is usual 
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to take the tabular points beginning at the point closest to 2* and taking 
subsequent points in increasing distance from x* so that 2* is near the 
middle of the values used. 


Example 2.4.1 (Computational) 

Use Neville’s algorithm to find an estimate for y(0.51) in table 2.1.1 
on page 9. Load INTERP and cycle on the Select method prompt 
to select Neville’s Algorithm. We shall use the option Direct so press 
RETURN again, and then enter the value z = 0.51. The values in 
table 2.1.1 are set as default values and should appear in a window on 
the right of the screen. We can now select the entries from this table in 
order of closeness to the point x = 0.51. The tabular value x = 0.5 is 
the closest one to z = 0.51, so that is the first one we wish to select. To 
do so, use the up/down cursor keys to position the inverse video block 
over the required value and then press RETURN. Another inverse video 
block will now appear under the heading ‘x, y selected’, and you can 
copy your chosen values into this just by pressing RETURN. The next 
tabular values needed are those at x = 0.7, which should be selected and 
copied in asimilar way. At this point type Q to stop selecting values and 
the program will then display a tableau based on the two values z = 0.5 
and x = 0.7. Further points can be added by typing Q and selecting 
the Select/edit data option. Add further values until you obtain the 
tableau shown in table 2.4.1. Notice that some editing of the selected 
values is allowed (should it be needed). For example, you can delete an 
entry (the entries below then move up one place), or copy one of the data 
points on top of another already selected entry which then replaces that 
entry. (Only three columns of the table of cross means are visible on the 
screen at any one time; other columns can be viewed using the left/right 
cursor keys.) Notice that by taking the tabulation points in such a way 
that the values |x* — 2;| are in increasing order of size the top entry in 
each column gives the interpolated value required. When the entries in 
successive columns are sufficiently close we can terminate the process. 
Additional values can be added to the bottom of the table in order to 
improve the current interpolated value. In this example the second and 
third entries agree to three decimal places and so it is sensible to quote 
the result as y(0.51) + 0.598. (In fact the value of po3 is much better 
than this; from section 2.1, P4(0.51) = 0.59734.) 


Exercise 2.4.2 (Computational) Use the program INTERP with the 
option Neville’s algorithm to estimate y(0.51) by adding the final entry 
to table 2.4.1 from table 2.1.1. 
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Table 2.4.1 Neville’s algorithm applied to table 2.1.1. 


z; (0.51 —2;) Yi Linear Quadratic Cubic 
0.5 +0.01 0.58813 


Po1 = 0.59733 
0.7 —0.19 0.77210 poz =0.59751 

pi2 =0.59367 pos = 0.59795 
0.3 +0.21 0.39646 pi3 = 0.6149 

p23 = 0.57133 


0.9 —0.35 0.89608 


Unfortunately, Neville’s algorithm does not always work as well as 
shown in example 2.4.1. If the tabulated data points are sufficiently 
smooth then all is well but consider the following example. 


Example 2.4.2 (Computational) 

Table 2.4.2 illustrates the application of the program INTERP, with 
option Neville’s algorithm, to estimate y(0.15), from the values given. 
The result is poor because the tabulated values are obtained from the 
function y(z) = 1/x which changes very rapidly near z = 0. (The 
correct value of y(0.15) = 6.666666, illustrating the difficulties which 
may appear if the data points vary greatly.) 


Table 2.4.2 Results produced by Neville’s algorithm. 


f y Linear Quadratic Cubic Quartic 
0.1 10.000000 
7.500000 
0.2 5.000000 7.083333 
5.833333 6.927083 
83% 3'333305 6.145833 6.848958 
4.583333 6.302083 
0.4 2.500000 5.208333 
3.750000 


0.5 2.000000 


Exercise 2.4.3 (Computational) Use the program INTERP to retrieve 
data from file SIN and use Neville’s algorithm to estimate y(47). The 
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data in this file is obtained from the function y = sinz. Insert addi- 
tional points and investigate how the accuracy of the interpolated point 
depends on the spacing and number of points. 


> 2.4.1 Inverse interpolation 


Thus far we have concentrated upon the problem of interpolating non- 
tabular values. In this section we will look briefly at the problem of 
inverse interpolation which was mentioned on page 10. The simplest 
solution to this problem is to reverse the roles played by the dependent 
and independent variables; for example, we simply rewrite table 2.1.1 on 
page 9 as table 2.4.3. This technique must be used with care because if 
there are two x values with the same y value then we will not get a single- 
valued interpolating function. This is equivalent to obtaining a singular 
matrix for equation (2.1.1). Furthermore, we need the original tabular 
values to be monotonic otherwise the order of the points will change 
when the roles of dependent and independent variables are interchanged. 


Table 2.4.3 Interchanging variables in table 2.1.1. 


Yield, y Concentration of reagent, z 
0.09956 0.10 
0.39646 0.30 
0.58813 0.50 
0.77210 0.70 
0.89608 0.90 


Exercise 2.4.4 (Computational) Use interpolation applied to table 2.1.1 
on page 9 to determine 2(0.4) using Neville’s algorithm and the op- 
tion Inverse. (Notice that applying inverse interpolation to this table 
is equivalent to applying interpolation to table 2.4.3. The program IN- 
TERP retains the original order of each (x,y) pair on the screen al- 
though their roles are reversed.) Check the result by using the option 
Construct interpolating polyn. and evaluating it at the point you have 
determined. How would you account for any discrepancy? 


The procedure of interchanging the dependent and independent vari- 
able can be very successful but must be used with care. It is clear that 
interpolation works best when the z; values are not too close and the y; 
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values do not change too quickly. If this is the case then we can expect 
a reasonably accurate interpolation, but if these conditions are not sat- 
isfied then inverse interpolation may be difficult. However, in order to 
illustrate the usefulness of inverse interpolation consider the following 
simple problem. 


Table 2.4.4 Inverse interpolation applied to x? — 3r+1=0. 


x y 

0.4 —0.040 
0.38261 

0.3 0.190 0.38189 
0.38636 0.38195 

0.5 —0.250 0.38261 0.38196 
0.39130 0.38213 0.38196 

0.2 0.440 0.38392 0.38202 
0.36296 0.38258 

0.1 0.710 0.37721 
0.34483 

0.0 1.000 


Example 2.4.3 (Computational) 
Use inverse interpolation to determine the solutions of the quadratic 
equation 

a — 321 =. 


Tabulating this function between z = 0 and x = 4 shows that the 
solutions lie in the intervals (0,0.5) and (2.5,3). Now load the program 
INTERP and enter values for f(z) = x? — 3x +1 for values of x 
between 0 and 0.5. Using Neville’s algortihm and the option Inverse 
produces the results in table 2.4.4. Notice that even though the roles 
of the dependent variable z and the independent variable y have been 
interchanged they are displayed in the same order on the screen. From 
this table the required solution is at c = 0.38196. (The solutions of 
the given quadratic are at c = $(3+ V5), the smallest of which is 
0.381966011.) 


Erercise 2.4.5 (Computational) Find the larger root of the quadratic 
equation given in example 2.4.3 by tabulating the function f(z) in the 
interval (2.5,3) and then using inverse interpolation. 
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Exercise 2.4.6 (Computational) Load the data from file SIN and use 
inverse interpolation to determine a value of z such that y = 0.5. The 
data in this file are obtained from the function y = sin 2; investigate the 
effect of additional points and different orderings. 


> 2.5 Spline interpolation 


Example 2.1.2 showed how it was possible to fit a polynomial of degree 
4 through five points, but the resulting function was quite oscillatory. 
We will now investigate how it is possible to fit lower order interpolants 
which are globally ‘smoother’. For example, through each pair of points 
(x;,y;) and (241, y¥:41) there is a unique straight line, but such a func- 
tion is unsatisfactory since it will in general have kinks at each data 
point. (See figure 2.5.1.) Mathematically speaking, the interpolating 
function will not be differentiable at the points z;,7=1,...,n—1. As 
an alternative, we could fit a quadratic function through (z;,y;) and 
(%i41, ¥i41)- This problem does not have a unique solution, but we can 
take advantage of this non-uniqueness to select the quadratic curves 
whose derivatives match at each of the data points. 


Figure 2.5.1 Local linear interpolation. 


Let us suppose that we attempt to fit the quadratic function 
Qi(x) =a; +b2+c27, fat ee Tt (2.5.1) 


in the interval (x;~1, 2;); then as there are n intervals we shall have 3n 
coefficients to find. First of all, we require that the resulting functions 
interpolate the points (z;,y;), 7 =0,...,n, Le. 


Q;(zi-1) = yi-1 and Q,(2;) = yi, tlie ny (2.5.2) 
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which gives 2n equations. (This condition will ensure that the resulting 
functions are continuous at the internal points z;,...,2n—1.) Further- 
more, at each internal point, z;, we require that the derivatives of Q: 
and Qi41 agree, 1.e. 


Oe y= Oa), fol eit ale (2.5.3) 
Therefore, we have 3n — 1 equations to determine the 3n coefficients aj, 
b;, cz, i= 1,...,n. It would appear that we can select any one coefficient 


arbitrarily. Combining together the resulting quadratic functions gives 
a quadratic spline. 


Figure 2.5.2 A quadratic spline through the data in table 2.1.1. 


The amount of effort in determining a quadratic spline can be re- 
duced by making use of the following procedure. Let us assume that in 
the interval (z;_1,2;) the required spline component is Q;(z), as shown 
in figure 2.5.3. We now assume that at z;_; the slopes of both Q;-1(z) 
and Q;(x) are d;_; and similarly at z;, Qi(z) = Qj,,(2) = d;. Since 
the derivative of a quadratic function is linear, Q/(z) is the linear cross 


mean of d; and d;_1, i.e. we have that 
(2 — aj21)d; = (2 a s)dji-4 
Li — Xji-1 


Q(z) = 


therefore, Qj(2i) = Q}41(2i) = di, which will ensure that the conditions 
(2.5.3) are satisfied. Next we integrate this expression to get 
_ (w@— aj-1)?d; — (x — xj)? d;_1 


Qi(z) = =i eo (eee Sy eee 
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Xi-4 Xj 


Figure 2.5.3 A quadratic spline component Q;(z) in (2;_1, 2;). 


where c; is a constant produced by the integration. Finally, we select the 
constants c; to ensure that the combined spline is continuous across each 
mesh node z = z;,2 = 1,...,n—1. (This ensures that equations (2.5.2) 
are satisfied.) Notice that 


OR y= ay Se RE oP 


and 
OA rpoa) 0s (2a 2s) dies + Cr = 4-1. 


We then eliminate c; between these equations to give 


2(yi-1 — Yi) 
d; = ——————~. -- d;_. 2.0.4 
Meva(@ici— 32) : Ooo 
This is a first-order recurrence relation and so all the required values, 
d,,dy,... can be found once do is specified. However, in general we 
will not be able to specify both dp and d,. (See Harding and Quinney 
(1986).) 


Example 2.5.1 (Computational) 

Load INTERP and when the prompt shows Select method press RE- 
TURN. Now change the option to Method: Quadratic spline and then 
press RETURN again. Use the graph option to set up suitable plotting 
ranges and then you will be requested to supply a value for the slope of 
the spline at its left end, i.e. the value do. For the moment accept the 
default value dy = 1.4845 and then observe the corresponding quadratic 
spline being plotted. (See figure 2.5.4.) Alternatively, we could take 
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dy = 0 which gives the quadratic spline shown in figure 2.5.2. The 
value do = (y1 — yo)/(z1 — Zo) is obtained by using a finite difference 
approximation for the slope at zo which explains why a better result is 
produced. (This is the default value which the program adopts.) 


Exercise 2.5.1 (Computational) Use the program INTERP to determine 
a quadratic spline through the data points in table 2.1.1. Investigate the 
effect of changing the value of do. (The quadratic spline for the default 
choice is given by figure 2.5.4.) 


Figure 2.5.4 Quadratic spline for table 2.1.1 with do = 1.4845. 


From figure 2.5.4 we see that a quadratic spline and its derivative are 
continuous at each node. However, the curvature is not smooth since 


: " d,—d 
ae Qo (2) i o = “ 
but 2 
Ate Qi(z Se men — _ 


and these may not be the same. A further problem with quadratic splines 
is that they tend to produce oscillatory functions as is also demonstrated 
in figure 2.5.2. In order to ensure that both the first and second deriva- 
tives are continuous we need to use a higher order local interpolant in 
each sub-interval. One of the most popular is to use cubic functions 
in each sub-interval and then attempt to select the coefficients so that 
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yj 
Mist 
Mj 
Xj Xia 


Figure 2.5.5 Cubic polynomial through (2;, y;) and (441, yi41)- 


both the slope and the curvature are continuous across internal tabular 
points. 

Let us suppose that we are given a table of data in the form of data 
points (z;, y;), 7 = 0,...,n. In the interval (z;,2;41) we use the cubic 
polynomial 

S;(z) = a; + bz + G2? + d,;x°. (2.5.5) 


In each sub-interval there are four coefficients to be determined, and 
since there are n such intervals this will give a total of 4n unknown 
values. Clearly, 


Si) ee and wee ert) age! — 10, ee 
which gives 2n conditions. Similarly, at each internal point we require 
Orn te) = Si(2z4), (eS il. 5 y= ins 


and 
Or = altel =a... al 


which gives a further 2n—2 conditions (see figure 2.5.5). It would appear 
that we can specify any two values arbitrarily. Once all 4n parameters 
are determined we shail have a unique cubic in each subinterval. This 
collection of functions is called a cubic spline. Even for a moderate 
number of tabulation points this procedure will result in a large system 
of linear equations but, fortunately, the total amount of effort can be 
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considerably reduced in a manner analogous to that used for quadratic 
splines. 
Notice that in the interval (z;, 2:41) we have that 


Si'(z) = 2c; + 6d;z. (2.5.6) 


Let us assume that the second derivative of 5; at 2; is given by some 
value M;, then 


Sa(2:) = 2c; + 6d;2; = 1k. 
CNS y) SEA Ey as eine 


Solving these equations for c; and d; gives 


Aes IVE eae VAs De Maina — Mia 2; 
: Lin. —2j : rhe ant 
Therefore, 
S!"(2) = (@ — i) Misi — (2 — tit) Mi 


Ti41 — Lj 


(This is just a linear cross mean.) Integrating this expression twice gives 


M; —z;)3 — M; — 2; 3 
S:(2) = Mele SAG a + One pe (2.5.7) 
where h; = 234; — z;. Now we fit this expression to the data points 
(a isy ian: (eres 1) ae: 
Mi(2i — 2:41)? + 6aiz; +68; = byi, 
Mi4i(ti41 — 21)? + 60i2i41 +68; = 6yi41, 
from which we eliminate 2; to produce 
oy = Ave) eee ie ;— M;) (2.5.8) 
(ti41 oH z;) 6 ( a+ 2 oO. 
and 
Bi = Yi — 25 — aMi(2i41 —2;)?. (2.5.9) 


However, S;(z) is required to be continuous across 2;, therefore 


j-1(21) = 5; (a), 
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M; (2; — 23-1) + 2aj-1 = Mi(2i — 2441) + 20; 
or 
Q@;—-Q;-1= $M; (2i41 — 2;-1). (2.5.10) 
Finally, we use equation (2.5.8) to eliminate a; and a;_; from (2.5.10) 
to give 


hy-1 Mj-1 + 2(hi-a + hi) Mi + AiMiai = S(it1 = Yi) _ B(yi = vi-a) 

h; Ai-1 
(2.5.11) 
This second-order recurrence relation holds at all internal points, there- 
fore, we have a system of n — 1 linear equations in the n + 1 unknowns 
Mo, M,,...,M,. In particular this system of equations has a special 
form since in a matrix notation we obtain entries only on the main di- 
agonal and the diagonals directly above and below this diagonal. These 
equations are therefore called tridiagonal and special methods have been 
developed to solve them which take advantage of their structure. (For 
an introduction to second-order recurrence relations and the solution of 
tridiagonal systems of linear equations see Harding and Quinney (1986).) 
The general solution of this recurrence relation will involve two arbitrary 
values which can be determined if two additional conditions are supplied. 
There are two particularly popular ways they can be specified: (i) free 
splines and (ii) clamped splines. 


(i) Free splines (natural splines). The simplest way of imposing ad- 
ditional conditions is to set both Mp and M, equal to zero. (This is 
equivalent to assuming that the interpolating function is approximately 
linear at zo and z,, since M; is the second derivative at z;.) 


Example 2.5.2 

Use the program INTERP to determine the natural spline through the 
data points given in table 2.1.1. As there are five points there will be 
five second derivative values to be determined. Let these values be Mo, 
M,, M2, M3 and Mg, but because the spline is to be natural we have 
Mo = Mz, = 0. Next we write down the equation (2.5.11) at the internal 
points, i.e. 


lo 


0.2Mo+0.8Mi+0.2M2 = (0.58813 —0.39646) _ (0.39646 —0.09956) 


j=) 
Li) 


|o 


© (9.77210—0.58813) — 
0.2 


6 6 
S(t) —0.77210) — —(0.77210 —0.58813 
Phe 89608 —0.77210) a ) 


0.2M,+0.8M2+0.2M3; = (0.58813 —0.39646) 


= 
to 


0.2M2+0.8M3+0.2M, 
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0.890: 25820 M, —3.156899 
0:20:89 072 My | = { —0.231000 |. 
OF 02255055 Mz —1.799699 


These equations have the solution M, = —4.306179, Mz = 1.440214 
and M3 = —2.609679. These values provide all that is necessary to 
construct the required natural spline. The resulting function is shown 
in figure 2.5.6. If we wish to evaluate the spline at a particular point it 
is first necessary to decide in which sub-interval the required point lies. 
For example, to evaluate this spline at x = 0.51 we need to evaluate 
the component of the spline which is defined on the interval (0.5,0.7) or 
(x2, 23). The required function is therefore 


or 


M3(z — x2) — Mo(z — 23)? 
= OS + fo, 
So(x) ee + apr + fo 
where 
a = al ll) SD bead 0.2(—2.660967 — 1.440214) 
0.2 6 
and 


fps 
Bo = 0.58813 — 0.5a2 — 5 * 1.440214(0.2)? 


from (2.5.8) and (2.5.9). Evaluating these expressions gives the compo- 
nent of the spline in the interval (0.5,0.7) and hence $2(0.51) = 0.59731 
to five decimal places. 


Exercise 2.5.2 (Analytical) Determine an explicit expression for the nat- 
ural cubic spline through the data points of table 2.1.1. By direct eval- 
uation show that the resulting spline is continuous and has continuous 
first and second derivatives at the internal points 2),...,2,_1. 


Exercise 2.5.3 (Computational) Use the program INTERP with the 
option Method: Cubic spline - free to confirm the results of example 
Pha? & 


(11) Clamped splines. If it is possible to supply values for the slope of the 
interpolating function at zo and z, then we can use equation (2.5.7) to 
provide the additional equations required. Let us suppose that we are 
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1.5 


1.0 


1.5 


Figure 2.5.6 The free cubic spline for the points in table 2.1.1. 


given that S¢(xo) = do and S/(z,) = d, then differentiating (2.5.7) and 
substituting we have 


- h 
gad iia Be G (Mi + 2Mo) 


ho 
or 
6 = 
2hoMy + hoM, = ote = Yo} =6do, (2.5.12) 
0 
and 
nn” Gn— hn- 
d, = (Uo —teo) + (Mk + 2Mn) 
n—-1 6 
or 
An—1My—1 + 2hn_1Mn = 6d, — gia Stoo) (2.5.13) 
n-1 


Next we combine equations (2.5.12) and (2.5.13) with the linear equa- 
tions given by (2.5.11) to obtain n + 1 linear equations in the unknown 
values M;, i = 0,...,n. Once these values are known the required cubic 
spline is uniquely determined. The matrix system which is formed in this 
way is both tridiagonal and diagonally dominant and can be solved in a 
relatively straightforward way by Gaussian elimination without partial 
pivoting. (See Harding and Quinney (1986).) 


Example 2.5.3 (Computational) 
Use the program INTERP to determine the clamped cubic spline which 
passes through the data points in table 2.1.1 and has slope 1.4845 at 
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z = 0.1 and 0.6199 at z = 0.9. (These values are default settings.) This 
clamped cubic spline is shown in figure 2.5.7 and can be compared with 
the free spline in figure 2.5.6. 


Figure 2.5.7 The clamped cubic spline in example 2.5.3. 


Exercise 2.5.4 (Computational) Use the program INTERP to deter- 
mine the clamped cubic spline which passes through the points given in 
table 2.1.1 and has zero slope at both ends. Compare the result with the 
free spline and the interpolating polynomial through the same points. 


We can now compare the quadratic spline which passes through the 
points in table 2.1.1, on page 9, which is shown in figure 2.5.4, the 
clamped cubic spline shown in figure 2.5.7 and the corresponding free 
cubic spline which is shown in figure 2.5.6. Clearly, the cubic spline is 
less prone to oscillations and we can show that this is generally the case. 


Theorem 2.1. Let C2[z0, tn] be the set of all twice differentiable func- 
tions which pass through the tabulation points, (z;,y;), 7 = 0,...,n. 
The natural cubic spline, S(x), which is a member of C2[z09,2,], mini- 
mizes the average curvature: 


i= : 7 if “ney az). 


Proof. Let S(x) be the required natural spline given by the solution 
of (2.5.11) subject to Mp = M, = 0 and let f(z) be any other twice 
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differentiable function over [xo, z,,]. Now consider 


/ “(rey de 


Tn 
ie) 


ii [S"(0) + (f"(2) — S"(a))}? az 


/ ” 9"(2)2 dx +2 i] ”” S(a)Lf"(2) — $"(a)]? de 


+ [U"e)- step aa 
and then look at the second of these terms, i.e. 
r= i S"f" — $") de. 
Integrating by parts gives 


i= S"(an)[f"(tn) = 9 (tn)] = S{2o)[f (20). —.9 (20)! 
— [ire -s(os"a) ae. 


However, for a natural spline S’”(ro9) = Mo = 0 and S”"(z,) = My, = 0 
and therefore the first two terms vanish. Furthermore, in each sub- 
interval (z;,2;41) the cubic spline has constant third derivative which 
we shall call K; and so 


t=n-1 


t=n-1 pri4gi 
a is [f'(2) - S'(a)| Kidz =— > KiLf(#s) - S(zi)] = 0. 


10) 


From this we conclude that 
[w@ra > [sera 


as required. 


In order to demonstrate the importance of this result consider the 
following example. 


Exercise 2.5.4 (Computational) Determine a polynomial approximation 
for the data in table 2.5.1. (These data points are stored in the file 
SPLINE). Now construct the natural cubic spline which passes through 
these points and compare the results obtained. 
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Table 2.5.1 Data for exercise 2.5.4. 


zx -—1 -—0.960 -—0.860 -—0.690 0.220 0.500 0.930 
y —-1 -—0.151 0.594 0.986 0.895 0.500 —0.306 


Erercise 2.5.5 (Computational) Use the program INTERP to deter- 
mine a quadratic spline through the data in the data file SPLINE which 
has zero slope at its left-hand end, i.e. dj) = 0. Vary do and investigate 
the behaviour of the resulting spline. Now fit a natural cubic spline 
through these data points and compare the results. Which curve is the 
smoothest? Investigate the behaviour of clamped splines for a variety 
of endpoint conditions. Is it possible to improve upon the natural cubic 
spline? To see that using a higher order interpolating function does not 
improve the behaviour, try fitting a high order polynomial using the 
option Construct interpolating polyn. (Hint: the quadratic spline is 
determined by a first order recurrence relation for which an initial con- 
dition is supplied. Is it stable? Why is the two-point boundary value 
problem formulation involved with cubic splines a more stable process?) 


> 2.6 Summary and conclusion 


We have seen in this chapter how it is possible to produce a continu- 
ous function which passes through a given set of tabular values. The 
solution to this problem is not unique. We can aim to fit a function 
which is composed of a linear combination of simple powers of 2, i.e. a 
polynomial, but this may lead to an interpolating polynomial which is 
highly oscilllatory. Alternatively, we can aim for a more local approxi- 
mation and fit a spline which will give a resulting function which is less 
likely to suffer from dramatic oscillations between tabular values. How- 
ever, whichever method is used we cannot say anything definite about 
non-tabular values. To illustrate this point see figure 2.6.1. The inter- 
polating function passes through all the tabular values but is likely to 
be of little use in interpolation. 


An underlying assumption with any interpolation problem is that the 
data points have been determined from some unknown but continuous 
function. If the underlying function is not continuous then no amount 
of continuous functions can produce satisfactory interpolation near any 
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Figure 2.6.1 An example of poor interpolation. 


singularity. Chapter 3 will consider this further, but for the moment we 
conclude this chapter with the following general remarks. 


1. Interpolation is usually more reliable than extrapolation; the latter 
should be avoided if at all possible. 


2. Except in all but the very simplest cases it is unwise to determine 
an interpolating polynomial by considering expressions of the form 
djp+a,;x+aozr?+...+a,2". This problem is possibly ill-conditioned 
and is likely to become more so as the number of tabular points 
increases. However, the construction of a polynomial interpolant 
can be carried out in a stable way using either Lagrangian poly- 
nomials or divided differences. If the value of the interpolating 
polynomial is only required at one or two points then it is better 
to use Neville’s algorithm. 


3. Interpolation is a process whereby given a table of values (2;, y(z;)) 
and a non-tabular value z* we attempt to find y(z*). Inverse 
interpolation is where we are given a value y and are asked to 
determine a value z which would produce y. 


4. The curve which has minimum curvature through a set of points 
and yet has continuous first and second derivatives is a natural 
cubic spline. Such a function is made up from a set of cubic func- 
tions, one in each tabular interval, which are continuous and have 
continuous first and second derivatives at each internal tabular 
point. 
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We have not considered the idea of rational interpolation, however, 
we shall consider the rational approximation functions in Chapter 3. 
Given a set of points and a rational function which passes through them 
it is clear that rational interpolation can be carried out at non-tabular 
points. 


> Chapter 3 


> The Approximation of Functions 


‘I’ve got fourteen pots of honey left, or is it fifteen, as the case may 
be’, said Pooh humbly. 
‘Well, let’s call it sixteen’ said Rabbit. 
With apologies to A A Milne. 


In an age of the ever increasing use of pocket calculators and micro- 
computers, the need to evaluate standard mathematical functions is also 
increasing. Clearly, electronic devices cannot store all possible values of 
such functions; the question therefore is how do such machines evaluate 
expressions such as sin(1.4) or exp(—3.456)? Inherent in such a problem 
is not only the need to carry out such an evaluation but also the need to 
know how accurate the presented result is going to be. In this chapter 
we will examine this question and try to provide some sensible answers. 


> 3.1 Taylor series approximation 


It is well known that provided a function is sufficiently smooth near a 
particular point x = a then it is possible to write 

z—a)? 

fle) = f(a) + (2a) f(a) + F= pra) +. 

where the infinite series converges for all z which satisfy |z—a| < Rand 
R is called the radius of covergence (see Appendix A). This is called a 
Taylor’s series expansion of the function f about the point z = a. An 
alternative form is 


es 2 
yay Gyn ea) f(a) p(ay tae 
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ae! = ar f(a) + Rn, (3.1.1) 


where the remainder term R,, is given by 


(x — ane aD) 
ee eG Be by 
and @ lies between z and a. Therefore, if we knew the value of f(a), 
f'(a), f(a), ..., f(a), we could use equation (3.1.1) to estimate f(z). 
Furthermore, the maximum possible error of such an approximation is 
given by max|R,| for all values of 0 between z and a. 


Example 3.1.1 

We can use a Taylor series with a = 0, which is sometimes called a 
Maclaurin series, to estimate exp(0.1). We can then consider questions 
such as the following. If five terms of this series are used what is the 
maximum possible error? How many terms are required to evaluate 
exp(0.5) correct to six decimal places? The Taylor series expansion of 
exp(x) about x = 0 is given by 


2 3 x4 n 


= : x Hb 
exp(z) = state “es pia sginsb ald taiteac eadia sean 


Taking the first five terms gives 


(Oa eal ae 
exp(0.1)%14+0.1+4+ 5 +: 6 8 4” 
therefore, exp(0.1) & 1.105170831 which can be compared with the cor- 
rect value of exp(0.1) = 1.10517092.... The remainder term R4 associ- 
ated with the above Taylor series approximation may be used to estimate 
the error involved with this approximate value. From equation (3.1.2) 
with z = 0 and a= 0.1 we have that 


0.15 exp(4) 
120 


Ral < e-exp(041) 
120 


m —38 
= o<reo.1 <9 2 Omer 


[sou 


compared with the actual error of 9.0 x 107. Suppose we now wanted 
to find exp(0.5) using the same quartic polynomial expression then 


0:57 0. 52a 025. 


exp(0.5) #14 0.5 + S-+ == + Fp = 1.6484374 
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but this compares less favourably with exp(0.5) = 1.6487213.... For this 
value of x, max |R4| = 4.29x 1074. Clearly, the range of values for which 
we can use this polynomial approximation is limited, but one way around 
this problem would appear to be to include further terms. Alternatively, 
we can use the remainder term, Ry, to see just how many terms are 
required. For example, in order to ensure that the approximation has a 
maximum error, for all |x| < 4, which is less than 1 part in 10° we need 
to select an n such that 


[Ral < 0.5 x 10-%, for all |z| < 2. 
Therefore, from (3.1.2) we need to find an integer n such that 


0.5°+1 exp(z) 


(a afl Se0 5 xe10 0 


|c|<3 


which can be achieved if n > 9. This result implies that if we require 
an accuracy of six decimal places each time we evaluate the exponential 
function for |z| < 5 then a polynomial of degree at least 9 is required. 
Clearly, for larger values of x more and more terms would be needed. 


Example 3.1.2 (Computational) 

We shall demonstrate a Taylor series approximation for the function 
exp(z), |z| < 1, using the program APPROX. (The program has a de- 
fault function setting for this problem.) Load the program APPROX 
and when the prompt Function is shown accept the default value as set 
which will produce a new prompt Method. Cycle through the options 
available by repeatedly pressing the SPACE BAR and when the option 
shows Taylor Series press RETURN. Again cycle through the options 
available, to familiarize yourself with the program, and when it shows 
Option: calculate coeffs press RETURN again. Now set the value of 
degree to 3, which will ask for a cubic approximation, and finally put 
a = 0. After a short pause the prompt will change to Plot: function; 
press RETURN to accept this. The prompt should now be Graph data: 
superimpose; however, clear the screen by changing this option to Graph 
data: clear — same axes and then press RETURN. The function will 
now be plotted in the graphics window and then the program will return 
to the option Plot: function. Now change the Plot option to approxi- 
mation by pressing the SPACE BAR and press RETURN to show the 
required cubic Taylor series approximation for exp(z) on this interval. 
It is also possible to plot the error between the function and approxi- 
mation by changing the plotting option, Plot, to show Plot: error. The 
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result is shown in figure 3.1.1. The associated error function is plotted 
on the same z axis but the corresponding y axis is scaled by a suitable 
factor which is displayed on the screen. In order to find this approxi- 
mation it is necessary to determine suitable values for the coefficients 
in equation (3.1.1) and with the option calculate coeffs this is done by 
using finite difference approximations. (See Appendix B.) Exercise 3.1.1 
demonstrates how it is possible for the user to supply these coefficients 
manually. 


Figure 3.1.1 Taylor series expansion for exp(z) about x = 0. 


Exercise 3.1.1 (Computational) Load the program APPROX and the 
function exp(z), |z| < 1. Select the Method: Taylor Series and the 
Option: calculate coeffs and then plot the function, its cubic Taylor 
series approximation and the associated error as described in example 
3.1.2. Exit from this part of the program by taking the option Plot: 
quit and follow the instructions given. This will return the program to 
the prompt Option: calculate coeffs. Again cycle through the options 
available at this point which include the possibility of storing the current 
coefficients on a suitable disc file, reading another set of coefficients 
which you have previously stored, or set coefficients manually. The 
degree of the current approximation should be 3 and the point about 
which the expansion is determined should be a = 0; if not reset these 
values. The correct values for the first four coefficients in the Taylor 
series of exp(xz) about z = 0 are 1, 1, ; and é. Accept the option to 
set coefficients manually and then input these values; notice that the 
current values, which are determined by finite differences, are shown 
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on the screen and are very close to the correct values. Compare the 
approximations which are determined. 


Erercise 3.1.2 (Computational) Use the program APPROX to deter- 
mine a cubic Taylor series expansion about z = 0 for exp(z), |z| < 2. 
Compare the results with a corresponding expansion about the points 
ee bean: ft, 


Ezercise 3.1.3 (Analytical) How many terms are required to evaluate the 
function log(1+ 2), |z| < 5, correct to eight decimal places? 


As the evaluation of polynomials plays a fundamental role in this 
type of approximation it is appropriate to consider the most efficient way 
of carrying out this task. We shall illustrate a technique called nested 
multiplication. This method permits the evaluation of a polynomial 
expression without the storage of intermediate results and therefore it 
can be carried out on a hand calculator which does not have a memory. 


Exercise 3.1.4 (Analytical) Polynomial evaluation by nested multiplica- 
tion. Let pn(x) = anz™ + G@,n2i2 +. + a1z + a9. Determine the 
number of multiplications and additions required to evaluate this poly- 
nomial directly for a fixed value of x. Show that if we define a sequence 
u; by 


Cp, Ts ede Se PST eer ees (3.123) 


then p,(x) = uo. This is called nested multiplication. Determine the 
number of operations required to evaluate the polynomial p, using this 
method. Which method is the most efficient? Notice that by using 
nested multiplication p,(z) can be calculated without any intermediate 
memory. 


For a fixed value of z the technique of nested multiplication intro- 
duced in exercise 3.1.4 is equivalent to a first-order recurrence relation. 
However, the recurrence relation (3.1.3) is stable only if |z| < 1 (see 
Harding and Quinney (1986)). To demonstrate the importance of this 
comment let us suppose that we wish to approximate sin(x), for a value 
of x in the interval [10,11], by using a Taylor series expansion of degree 
3 about the point z = 10. This gives 


p3(z)=sin(10) + (x—10) cos(10) — 4(x—10)? sin(10) — 4(x—10)° cos(10). 
(3.1.4) 


Expanding out this expression gives the cubic 


p3(x) = ao + ayx + aor? + a32° for x € [10,11], (3.15) 
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where ay = —104.798, a; = 35.6743, ag = —3.92335 and ag = 0.139845, 
correct to six significant figures. Now suppose we wish to find p3(10.5). 
As a3 is accurate to six decimal places, when it is be multiplied by a 
factor (10.5)? the resulting value of p3(10.5) may only have two or pos- 
sibly three decimal places correct. Clearly, the evaluation of the poly- 
nomial (3.1.5) may be prone to induced errors and should be avoided. 
Instead we leave it in the form (3.1.4), i.e. 


p3(x) = —0.544021 — 0.839072(¢ — 10) + 0.272011(z — 10)? 
+ 0.139845(x — 10)° 


which ensures that all the terms raised to powers are less than 1; hence 
we should have a stable recurrence relation and reasonable condition- 
ing. This introduces a very important point, for clearly this problem 
will always arise if we attempt to approximate any function outside the 
interval either [0,1] or [—1,1]. However, we can always scale any inter- 
val [a,b] to either of these standard intervals by a suitable linear scaling 
procedure. 


Exercise 3.1.5 (Calculator) By examining the possible errors in the coef- 
ficients of the polynomial (3.1.5) determine upper and lower error func- 
tions for this expression. 


Exercise 3.1.6 (Analytical) Scale the interval [2,10] to the standard in- 
tervals [0,1] and [—1, 1]. 


> 3.2 Polynomial approximation 


We have seen how it is possible to obtain a polynomial approximation 
for a given function by truncating a corresponding Taylor series expan- 
sion. Now we will consider an alternative approach using interpolating 
polynomials. Let us suppose that we are given a function f(x), z in 
(a,b), and attempt to find a polynomial function of degree n, Dal ts 
which is close to f in some sense. 


Example 3.2.1 (Computational/Calculator) 

Let us determine a polynomial approximation of degree 3 for the function 
f(x) = exp(z), for values of z in [—1, 1]. We begin by considering the 
Taylor series expansion of exp(z) in this interval. This is given by 


exp(z) =14+2+4 $274 zx° + R3(z), 
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where the remainder term satisfies |3(xr)| < maxyz)<1 |f*?|/4! < 0.113. 
This approximation is shown in figure 3.2.1 and the error function, 
Ra(z), in figure 3.2.2. As an alternative method we can make the func- 
tion and approximating polynomial agree at a given set of points, i.e. we 
select the coefficients a nth degree polynomial so that pn(z;) = f(2;), 
27=0,..,n+1. Thus the problem becomes one of constructing an in- 
terpolating polynomial through n+ 1 points. Firstly, we need to decide 
where we are going to put the tabulation points z;, 1 = 0,1,2,3, at which 
the polynomial and function agree. For the moment, let us put them at 
equally spaced points in [—1, 1], i.e. —1, 3,3 and 1, and assume that 
the required cubic polynomial is 


p3(x) =agtayxr+ ayx? + a3x°. 


Then at 2 = —1 we require that 
p3(—1) = do\— 4; + ay — ag = exp(—1). 


We obtain a similar equation at the remaining three points which gives 
a system of four linear equations in the unknown coefficients ag, a1, a2 
and az (see section 2.1). These equations can be solved using Gaussian 
elimination to produce the cubic approximation 


p3(x) = 0.995196 + 0.999049x + 0.5478852? +. 0.1761522°. (3.2.1) 


This function is also plotted in figure 3.2.1 and its error in figure 3.2.2. 
Notice that the maximum error associated with the polynomial (3.2.1) 
is much less than the maximum error associated with the corresponding 
Taylor expansion of the same order. 


We have already seen how it is possible to estimate the error asso- 
ciated with a Taylor series approximation of a given function. Now we 
will construct an error function for polynomial interpolation. Firstly, we 
recall that f(x) and the interpolant p,(z) agree at the n + 1 points 2;, 
(— lea, ity bnerelore, 


én(z) = f(z) — pn(z) = K(e)(x — z0)(z — 21) ...(2 — Zn), 


where K(z) is to be determined. For any given value of z let us define 
the function ¢ by 


P(t) = f(t) — prt) — K(z) w(t), 
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Figure 3.2.1 Taylor series approximation and p3(z) for exp(z). 


4 Error in cubic 
Taylor series 


Error for cubic 


Figure 3.2.2 Error functions associated with figure 3.2.1. 
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where w(t) = (t—zo)(t—21)...(t—2,), and then notice that ¢ vanishes 
at the points n + 2 points x, x9, ...,2n. By Rolle’s theorem ¢’ must 
vanish between any two points at which ¢ vanishes and so ¢’ has n+ 1 
zeros. Similarly, $” has n zeros and so on until we find that the (n+1)th 
derivative of ¢, ¢("+), has exactly one zero; let this zero be 6, then 


0 = gOrFN(B) = FON (6) — pr+D(A) — Kw (6), 


However, prt) is identically zero and w\"+!) = (n + 1)!, therefore, 
(n+1) 
teres A!) Ge 
(n+1)! 
or 
i fA) 25% 
PITA? Wess 19 |G 30 a Ce are ( 4G ) 


(Recall that 6 will depend upon the choice of z.) Notice the similarity 
between the form of this error function and that associated with the 
Taylor series expansion given by equation (3.1.2). 


Example 3.2.2 (Calculator/Computer) 
We shall determine an upper bound on the error associated with the 
cubic approximation for exp(x) which was derived in example 3.2.1, i.e. 


p3(z) = 0.995196 + 0.999049x + 0.54788527 + 0.1761522°. 
From equation (3.2.2) with n = 3 we have that 
f(x) — pa(2)| < |(e + 1)(@ + 3)(@ — 3)(@ — 1)[M/A!, 


where M = sup |f*’|, x € [—1,1], and f(z) = exp(x). (We define the 
supremum as follows: if M > |f*”| and V L > |f?**|, M < L then M is 
the supremum.) We now examine the first four terms and note that for 
values of z between —1 and 1 this takes a maximum value of 0.19753 
ah ed +15. (See figure 3.2.3.) Therefore, since |f*’(x)| = |e”| < e, 
for x € [—1, 1], the error function e3(z) satisfies |e3(x)| < 0.0224, which 
is the required error bound. In fact this is a severe overestimate since 
the maximum error is actually 0.00481 but even a crude bound is better 
than none. (Recall that the corresponding bound on the cubic Taylor 
series approximation was |R3(zr)| < 0.113.) 
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0.4 


0.3 


Figure 3.2.3 y = |(e+ (e+ H(e-H(e- Dh, lel <1. 


Exercise 3.2.1 (Analytical) By determining upper and lower bounds on 
the fourth derivative of exp(z) construct upper and lower bounds for 
the error function e3(z) given in example 3.2.1. 


Exercise 3.2.2 (Computational) Use the program APPROX to deter- 
mine a polynomial approximation of order 2 for exp(zx), |z| < 1, by inter- 
polation at the points z = —1, —3 and 1. Load the program APPROX, 
accept the option Function as set and then select Method: Interpola- 
tion. Now select Option: edit /inspect points and enter the three points 
on the curve at z = —1, —} and 1. (When you have selected this option 
you will be given the opportunity to change the plotting ranges if you 
wish to do so.) Select the editing option, cycle through the edit options 
until it shows Add Point and then press RETURN. The top window 
will now contain information regarding how to input a point and, in the 
graphics window, there should be a flashing cursor at the point (—1,0). 
For example, if you press the COPY key the cursor will move to the 
point (—1,exp(—1)). If you press RETURN this value will be accepted 
as an interpolation point and the number of points, N,, which is shown in 
the bottom window will be increased by 1. Finally, the edit option will 
be returned to Add Point. Now input two additional points at xz = -3, 
ie exp(—4) and z = 1, y= exp(1). The z values can either be input 
by pressing ‘X’ or by using the cursor which is controlled in the usual 
way. (Each key press will move the cursor by an amount ‘dX’ and this 
variable is in turn controlled by the ‘>’ and ‘<’ keys.) When you are 
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satisfied with the points which have been input quit the editor and select 
Option: calculate coeffs then press RETURN. The computer will now 
determine the polynomial through these points, display the resulting 
polynomial in the bottom window and then display the prompt Plot: 
function. By changing the prompt you can plot the approximation or 
the error between the function and the computed polynomial. Use the 
editor to change the interpolation points and investigate how the approx- 
imation and its error change. Now return to the editor and add another 
two points and then determine quartic polynomial approximations for 
exp(z), x € [—1, 1]. 


Erercise 3.2.4 (Computational) Use the program APPROX with the 
option Method: Interpolation to determine the cubic polynomial ap- 
proximations for the function exp(z), x in [—1,1], by interpolation at 
—1, —0.5, 0.5 and 1. Investigate how the error between the interpolat- 
ing polynomial and the function changes as the tabular points change. 
Investigate how the error changes as the number of tabular points in- 
creases. Find the polynomial of smallest degree which approximates this 
function to an accuracy of at worst 0.001 throughout this interval for 
interpolation points which are equally spaced. Investigate how the error 
function behaves for different interpolation points. 


We have assumed that it is always possible to construct a suitable 
polynomial approximation, but is this always the case? Fortunately, 
provided the function f satisfies certain minimal requirements then there 
is always such a polynomial. We state, without proof, the following 
theorem. 


Theorem 3.2.1. Weierstrass’s approrimation theorem. If f is a continu- 
ous function on the finite interval [a, 6], then for any € > 0 there exists 
an nth degree polynomial p,(x), for some n, such that 


|f(x) — pn(x)| < € for all x € [a, 6). 


In other words, if we are given a continuous function f, defined on the 
interval [a,6], then there is always a polynomial function which will 
approximate f to within any prescribed error €. Unfortunately, what 
Weierstrass’s theorem does not tell us is how large n should be, nor how 
to find the associated polynomial. Furthermore, this theorem assumes 
that we can carry out all arithmetic exactly, which is not the case either. 
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> 3.3. Minimax polynomials 


In the last section we saw that there are many polynomial approxima- 
tions of a given degree for any given function. For example, 


exp(z) ¥1+2+4 hy? + iz, (3:3.1) 
and 
exp(z) © 0.99520 + 0.999052 + 0.547882? +0.1761523, (3.3.2) 


are both cubic approximations. The first is a truncated Taylor series, 
the second is an interpolating polynomial which agrees with exp(z) at 
the points —1, —5, ; and 1. We can compare these functions as shown 
in figure 3.2.1 but it is more informative to iook at their respective error 
functions shown in figure 3.2.2. If we were to obtain another interpolat- 
ing cubic approximation by adjusting the tabular points we could plot 
the associated error function in the same way and compare it. In gen- 
eral we need to find a polynomial, p,(x), which minimizes the maximum 
value of the error function, i.e. for any other polynomial of degree n, 
Qn(z) say, the minimax error €,, satisfies 


€n = max|f(z) — pr(x)| < max|f(xr) — qn(zx)| for all x € [{a, 6]. 


The resulting polynomial is called the minimaz polynomial of degree 
n. By adjusting the tabular points at which a polynomial q,(z) agrees 
with f(z) we can look for the minimax polynomial, but this is likely 
to be a long process. Instead we can use a very powerful result due to 


Chebyschev. 


Theorem 3.3.1. Chebyshev’s equt-oscillation theorem. Let f(x) be a 
continuous function defined on the interval [a,b] and p,(x) be a poly- 
nomial approximation, of degree n, for f on [a,b] with associated error 
function en(z) = f(x) — pr(z). If there exists a set of at least n + 2 
points, a < 2% < 21 <... < Xn41 < b, where the error function takes 
maximum values such that 


Cn fy ee Ea... C10) 18. yr 


then pp(z) is the minimax approximation. This approximation is the 
only one of degree < n with this property. 


MINIMAX POLYNOMIALS 59 


Notice that the equi-oscillation theorem only requires alternating 
minima and maxima in the error function; they need not be turning 
points as we shall see in the next example. This theorem is extremely 
important for it enables us to decide whether or not a particular polyno- 
mial is minimax amongst all polynomials of a given degree. For example, 
the error associated with the third-order approximation (3.3.2), which 
is shown in figure 3.2.2, has five alternating minima and maxima at 
the points P, Q, R, S and T, but since they are not all of the same 
magnitude this polynomial is not minimax. 


= 0 1 


Figure 3.3.1 Linear minimax approximation of exp(z), |z| < 1. 


Example 3.3.1 (Analytical) 

In order to demonstrate the minimax equi-oscillation property we will 
construct a minimax linear approximation for exp(x) on the interval 
[—1,1]. Let us write the linear minimax approximation for exp(zx) as 
p(x) = ax + 5, then the error function e;(z) is 


€,(z) = exp(z) — az — b. 


As the degree of the polynomial approximation is 1 we need to determine 
three points x9, 2; and 2x2 such that e;(x) takes maximum absolute 
values which are of the same size but alternating sign, 1.e. 


€1(z0) = E, e1(21) = —E, and e;(z2) = E. 


At least one of x9, 2; and z2 must be an internal point of the interval 
[—1,1] at which there will be a turning point, ie. de,/dx = 0, which 
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gives 

exp(z) -a=0. 
This equation has only one solution, therefore only one of zo, 2; and £2 
is a turning point. Let this be z,. Accordingly, e;(z) will have maximum 
absolute value at both zo and z2 neither of which will be turning points. 
Therefore, 


exp(-l)+a-b = E 
exp(z1)-—az, -b = -E 
exp(l)-—a-b = E 


which gives three equations in the four unknowns; the unknown point 
21, the coefficients a and b, and the error FE. However, at x, the error 
function has a turning point and so 


which gives a fourth equation. The solution of these equations is: 


(exp(1) — exp(—1)) 


: = 1.17520119, 


zy) = loge(a) = 02161439361" 
b {exp(ei) + exp(}) ~ 971 — a) — 1.26427905 
E = exp(1) —a— b = 0.278801586. 


The resulting polynomial approximation is shown in figure 3.3.1 together 
with its error function in figure 3.3.2. Recall that we are seeking a 
polynomial of degree 1 and notice that the error function in figure 3.3.2 
has three alternating minima and maxima at P, Q and R, all of which are 
equal to FE in magnitude. Therefore, this function is the linear minimax 
approximation for exp(x), for z in [—1,1] and E is the minimax error 
€1. In order to investigate this idea further consider an alternative linear 
approximation for exp(z) given by q:(z) = 1+ 2, the error function for 
which is also shown in figure 3.3.2. It is clear that the maximum error 
of q; is greater than that of p;. The Chebyshev equi-oscillation theorem 
states that for any other linear approximation the associated error has 
a maximum value greater than the error associated with the minimax 
linear approximation p;(z). For linear minimax approximations it is 
possible to determine the correct expression algebraically, but in general 
this is not the case. 
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Figure 3.3.2 Linear minimax error function for exp(z), |z| < 1. 


Example 3.3.2 (Analytical/Computational) 

We shall attempt to construct a cubic minimax approximation for exp(z) 
for x in [—1, 1]. Firstly, let us examine the problems of determining the 
minimax approximation analytically. We need to find a set of five points, 
(n = 3 since we seek a cubic approximation), ro, 21, £2, 3 and x4, where 
the error function 


e3(z) = exp(z) — agz® — aor? — a,x — ag 
has alternating minima and maxima all of which are of equal magnitude. 
Differentiating this function shows that de3/dz has no more than three 
zeros in [—1,1] and so zg = —1 and z4 = 1 at which the error function 
will have maximum absolute value. (To see this, sketch the curve y = 
exp(z) and any other cubic.) Therefore, we have to solve the equations 


exp(—1)+a3-—a2+a,;-a9 = E 
exp(2,) — azz? — doz? —a,2,-a) = —E 
exp(Z2) — a3x3 — a2x> —1j2o—d9 =. E (3.3:3) 
exp(x3) — a3z3 - aor? —a,|23—a = —-E 
exp(1)-—a3-—ag-—a,;-a)9 = E 


plus the conditions for turning points at 21, 22 and 23, 1.e. 


exp(z1) — 3a3x7 —2a0z, —a; = 0 
0 


exp(x3) — 3a3%3 —2a9z3 a, = 0 


exp(x2) — 3a323 — 2agr2 — ay 
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which gives a system of eight non-linear equations in the unknown values 
21, £2, £3, the error EF, and the coefficients ao, a), a2 and a3. In general 
it would be extremely difficult to solve these equations. 


Exercise 3.3.1 (Analytical) Construct a cubic minimax approximation 
for the function cos(x) for values of x in [—1,1]. (Hint: use symmetry 
to reduce the size of the problem.) 


Example 3.3.2 demonstrates the difficulty in constructing minimax 
approximations of order n > 1. We will now consider an alternative 
iterative approach. Consider the cubic approximation for exp(z) given 
by (3:241)pi-e 


exp(xz) © 0.99520 + 0.999052 + 0.54788x” + 0.176152°, 


which is shown in figure 3.2.1. The associated error, as shown in fig- 
ure 3.2.2, has the correct number of minima and maxima but they are 
of different magnitude. The internal stationary points of the error func- 
tion occur at approximately —0.7, 0 and 0.8, therefore let us solve the 
equations (3.3.3) with 2; = —0.7, z2 = 0 and x3 = 0.8. This produces 
the linear equations 


1 -1 1 —1 1 dao exp(—1) 
1 -0.7 0.77 -0.77 —-1 ay exp(—0.7) 
Lise at) 0 0 1 oy | exp(0) 
14750, 85¢ 50:82: e087 ame a3 exp(0.8) 
1 1 1 1 1 E exp(1) 


and in turn the solution of these equations gives the cubic 
q3(z) = 0.99465 + 0.996432 + 0.543082? + 0.178772°. (3.3.4) 


The error e3(x) = exp(z) — q3(z) is shown in figure 3.3.3, together with 
the computed value for FE, and has turning points at —0.687312, 0.045814 
and 0.727528 in addition to the maxima at —1 and +1. (These extremum 
points are denoted by P, Q, R, S and T). As we seek a minimax 
approximation of degree 3 the error should have 3+2 alternating extrema 
of the same magnitude. From figure 3.3.3 we can see that the error 
function does have five extrema which have alternating sign but they are 
of different magnitude. If the extrema were of the same magnitude then 
E would be the minimax error €3 for exp(x). We could now substitute 
these five points into (3.3.3) to determine another approximation and 
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hence iteratively refine our cubic approximation. (Notice that only the 
zx co-ordinate is required to determine the next approximation.) This is 
the basis of the Remes algorithm for determining minimax polynomial 
approximations which can be described schematically as in figure 3.3.4. 


Figure 3.3.3 Error curve associated with the cubic (3.3.4). 


Example 3.3.3 (Computational) 

We demonstrate the Remes algorithm by considering a cubic minimax 
approximation for exp(x), —1 < z < 1. Load the program APPROX, 
for which this is the default function. When the prompt shows Function: 
as set press RETURN and select the option Method: Minimax/Remes. 
As we seek the minimax approximation of order 3, we will need to sup- 
ply five points at which the required approximation for exp(x) has al- 
ternating minima and maxima. For the moment let us use the values 
xz = —1, —0.7, 0, 0.8 and 1 which can be input using the prompt Option: 
edit /inspect points. When you are satisfied that suitable points have 
been entered quit this option and change the prompt to Option: calcu- 
late coeffs. (It is also possible to supply a first approximation for the 
minimax polynomial using the option Taylor Series). Now press RE- 
TURN and the computer will determine the solution of the equations 
in step 2 of the Remes flowchart, shown in figure 3.3.4, and display the 
prompt Plot: function. If the function has not been plotted prior to 
this point then pressing RETURN will display it in the graphics window, 
otherwise, change the prompt to Plot: approximation and compare the 
computed approximation with the function. It is likely that the function 
and the computed approximation are very close; the error between them 
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1. Select m+2 points 
BG Oy 6 oe 


iXpers 


2. Solve the linear equation 
p,(x,)-#(x; )=(-1)'E for the 
coefficients of p, and E. 


3. Determine the minimum and 
maximum points of the curve 
B(x) - F(x). 


. Are there 
exactly 7+2 
extrema? 


MES 


6. Select n+2 of the extrema 
such that the largest extrema 
are selected and that they have 
alternative sign. Label these 
PORTS Exon ee 


5. Label these points 
XQueis aX piace 


7. Test for 
convergence 


Figure 3.3.4 The Remes algorithm. 
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can be examined by using the option Plot: error. For this example you 
should obtain figure 3.3.3 together with a dotted line which denotes the 
value for EF in equations (3.3.3). The coefficients of the computed ap- 
proximation can be examined by using prompt Option: tabulate, which 
in this case gives the cubic (3.3.4). Notice that the error associated with 
this approximation has maximum values at the points P, Q, R, S and 
T’. We can select these points to be used in the next stage of the Remes 
algorithm by using the prompt Option: edit/inspect points and then 
return to Option: calculate coeffs to determine the next approxima- 
tion. (Only the x coordinates are required, therefore, it is unnecessary 
to enter the y coordinates.) This process can then be repeated until 
convergence is obtained. 


Ezercise 3.3.2 (Computational) Use the program APPROX to carry 
out the next iteration of the Remes algorithm using the minima/maxima 
of the error function associated with the cubic approximation (3.3.4). 


Ezercise 3.3.3 (Computational) Use the program APPROX with the 
option Remes/minimax to determine the quartic minimax approxima- 
tion for exp(x), |z| < 1. 


Exercise 3.3.4 (Computational) Use the program APPROX with the 
option Remes/minimax to determine the quartic minimax approxima- 
tion for log(1 + 2), |z| < rt Try different starting polynomial approxi- 
mations to examine the speed at which this algorithm converges to the 
required minimax approximation. 


At each stage of the Remes algorithm it is to be expected that the 
maximum value of the associated error function will be reduced but it 
will certainly be an upper bound on the minimax error €, for polynomial 
approximations of order n. In fact we can also produce a lower bound 
for the minimax error €y. 


Theorem 3.3.2. (de la Vallee Poussin). 

Given that g, is an nth degree polynomial approximation for a function 
f(z),x € [a,], and z;, i = 0,...,n +1, are any n+ 2 points in [a, }] 
such) thatia,< Zoom 2, < .-) < <n41 < 6 where the error function 
€n(x) = f(x) — qn(x) alternates in sign at successive points, then 


Tene Nea Cn aX. | (. 2); 
Ly xr€[a,b] 


where €,, is the minimax error. 
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Notice that the upper bound in theorem 3.3.2 is over the whole in- 
terval, but that the lower bound is only at the points xo, 41, ..-, n41- 
For example, the cubic approximation (3.3.4) produces the bounds 


0.0053526 < e, < 0.0059120. 


(The lower bound is taken over the points —1, —0.7, 0, 0.8 and 1.0 and 
occurs at z = —0.7. The upper bound is over the whole of the interval 
and occurs at x ~ 0.7275.) If these bounds are sufficiently close to each 
other then there is little point in carrying out another iteration of the 
Remes algorithm. What is more, the lower bound tells us that if we wish 
to have a minimax error less than 0.0053526 then we must construct a 
polynomial approximation of higher order. 

The Remes algorithm is based upon the selection of exactly n + 2 
alternating minima/maxima from a given error function, but it is pos- 
sible to obtain an error curve for a polynomial approximation of degree 
n which has more than n+ 2 such points. (Although p,(z) cannot have 
more than n — 1 turning points f(z) — p,(z) could have any number.) 
In this case we simply select n + 2 values where the error function has 
alternating minima and maxima, the remainder are discarded, the only 
proviso is that we must take the largest minima/maxima amongst the 
points we select. We will not give details of how to find the relevant 
minima/maxima but in general this is usually accomplished by a linear 
search along the error curve, not by differentiating. However, the loca- 
tion of suitable points by hand is a relatively simple procedure as the 
following example demonstrates. 


Example 3.3.4 (Analytical) 

Let us suppose that we have constructed a quadratic approximation for 
some function and the error curve looks like figure 3.3.5. The minimax 
approximation should have four alternating minima and maxima, this 
curve has five. We need to select four of the values shown in order to 
proceed with the Remes algorithm; in this case we take 21, 29, 23 and 
rq as these are alternating minima and maxima and the maximum error 
is at x4. The point 29 will no longer be used. 


As the order of the polynomial approximation increases the time 
required by the Remes algorithm increases rapidly and the closer the 
initial approximation is to the minimax the better. Furthermore, as 
the number of internal points increases equation (3.3.3) may become 
ill conditioned. The reason for this is clear. As n increases then the 
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Figure 3.3.5 A possible quadratic approximation error curve. 


matrix corresponding to equations (3.3.3) becomes very close to the 
Vandermonde matrix discussed in Chapter 2 which is known to be ill 
conditioned as the tabular points get closer together. A solution to this 
problem will be discussed in a later section. 


Exercise 3.3.5 (Computational) Use the program APPROX to con- 
struct a minimax polynomial approximation for exp(xz) which is of min- 
imal order but has error no worse that 0.001 for -l <2 <1. 


Exercise 3.3.6 (Computational) Use the points —1, 0, 1 to construct a 
linear minimax approximation for x° over the interval [—1, 1]. Explain 
why the Remes algorithm fails in this case and suggest how it is pos- 
sible to overcome this difficulty. (Hint: look at the quadratic minimax 
approximation.) 


> 3.4 Chebyshev approximation 


All the approximations we have discussed so far have been simple power 
series approximations, i.e. of the form 


i=n 
(EARS esa hewr ee seas (3.4.1) 
i=0 


in the hope that, provided the series converges as n — oo, taking a fi- 
nite number of terms will provide a suitable approximation. If we wish 
to obtain an approximation close to the minimax by truncating equa- 
tion (3.4.1) after n terms then the remainder must satisfy the Chebyshev 
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equi-oscillation property. Clearly, using this simple power series cannot 
produce this phenomenon, therefore, we will now look at alternative 
expansions. We begin by looking at Chebyshev polynomials. 

Let us consider a sequence of functions which are defined, for all 
values of x € [—1, 1], as follows: 


Tye) (2) —z, 
Tae (ya ed ae ele ey, (3.4.2) 


From the recurrence relation (3.4.2) we can generate 


2xT; (x) — To(z) = 2x? — 1, 
22To(x) — Ti(2) = 42° — 3z, 


T2(z) 
T3(z) 


and so on. These functions have some remarkable properties. For exam- 
ple, notice that 7;,(z) is a polynomial of degree n and that for n odd T, 
is an odd function whilst for n even T;,, is even. The functions 7;,,(x) are 
called Chebyshev polynomials. The first four are shown in figure 3.4.1. 


Figure 3.4.1 The Chebyshev polynomials 7;(z),i = 0,1,2,3. 


Ezercise 3.4.1 (Analytical) Show that T,(z),z € [1,1], is a polynomial 
of degree n and that it is odd when n is an odd integer and even when 
n is an even integer. (Hint: use induction.) 
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We can also define Chebyshev polynomials in another way, i.e. 
Te (2) cos( cosa (2). (3.4.3) 
To see this note that 


Tn4i(z) = cos((n +1) cos~}(z)) 
= cos(cos~'(x)) cos(n cos~(x)) 


—sin(cos~‘(xr)) sin(n cos~}(z)) 
and 


Tee cye =" cos(( —-1)*cos + (2)) 
= cos(cos~!(x)) cos(n cos +(x) 


+ sin(cos~!(z)) sin(n cos~(z)) 
Adding together these equations gives 


Tn41(2) + Tp—1(2) = 2 cos(cos~+(z)) cos(n cos~!(x)) = 2rT;(z), 
(3.4.4) 

which is equivalent to (3.4.2). The definition (3.4.3) now enables us 
to draw another conclusion about the polynomial 7;,(z): it has exactly 
n+ 1 alternating minima and maxima, with equal magnitude, in the 
interval [—1,1]. The importance of this comment can be demonstrated 
as follows. 

Let us assume that instead of (3.4.1) we can write 


=n 


i 
fay =e cl (2)-b eneiinst + entoInge +. <<. (3.4.5) 
i 


{| 
fo) 


If this is the case then the dominant term in the error when this series 
is truncated after n + 1 terms will be cn417n41 and as 7,41 has ex- 
actly n+2 equal and opposite minima and maxima we should obtain an 
approximation close to the minimax! We must now justify the expan- 
sion (3.4.5) of any function in terms of Chebyshev polynomials. First 
of all, if the series (3.4.1) converges for z € [—1,1] then so will any re- 
arrangement of its terms, but this is just what the Chebyshev expansion 
is. Therefore, if (3.4.1) converges then so will (3.4.5). The next stage is 
to find an expression for the coefficients c;, 7 = 0,1,.... This turns out 
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to be quite straightforward because the Chebyshev polynomials possess 
an orthogonality property, i.e. 


ee on 
SU rs aor hms 
Pen/ ( (1 = pe 0 
To see this notice that the substitution z = cos@ reduces 
ia) late) Fae 
1 Y(1- x2 


to the standard trigonometric integral 


“ eet — Sie () 
| cos 76 cos 78 dé = i=) 0 
0 0, i#;3 


Therefore, in order to determine the coefficients c; in equation (3.4.5) we 
multiply both sides by T;(x)/(1—2x?)!/? and integrate over the interval 


(—1, 1), i.e 


(3.4.6) 


ie resers =f Hleoe8) 8, 


a = eee cos cos 
aye vet Se) 5 (6)) cos j0 d0. (3.4.7) 


Alternatively, we can redefine the Chebyshev series (3.4.5) as 


t= Co 


= eo + y oTi(a)+ CrsiTnait Cn4oin4o tee ee 


ar 

’ (3.4.8) 
where the symbol denotes that the first coefficient is i: All the 
coefficients are then given by (3.4.7). 


Example 3.4.1 (Analytical) 
Determine the Chebyshev coefficients for the function 32(1 — 2?)!/?. 
Notice that this is an odd function, i.e. f(x) = —f(—z), and so all the 


even coefficients will vanish; we need only determine the odd coefficients. 
Therefore, 


cotees 2f 32,/(1 — 22)T(z) ate 


(bee 2d) 
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4 1 
= =| 327? dz = 
T Jo 


ee 2f 32/(1 = 2#) Ta(z) | 
LEA Es) Y(1 — x?) 
12 


4 1 
= =a 32(42° sae 32) dz = =) 
0 on 


Tv 


) 


= [oe 


and so on. The function given by ¢,7;(r) + c373(x) is compared with 
f(x) = 3x(1 — x?)1/? in figure 3.4.2 and the associated error is shown 
in figure 3.4.3. Notice that the error function has six alternating min- 
ima/maxima, but of unequal size. Although we have approximated f 
by a cubic, the six minima/maxima cause no problems because the co- 
efficient c4 is zero and therefore the quartic approximation for f is given 
by a cubic. 


1.0 


beer 
3xV(1—x*) 


Figure 3.4.2 Cubic Chebyshev approximation for 32(1 — ee 


Example 3.4.2 (Analytical) 

In order to illustrate approximation using Chebyshev polynomials we 
shall compare a linear Chebyshev approximation for the function e”, 
|z| < 1, with the corresponding linear minimax approximation Dyke) 
1.2464279 + 1.175201z, which was determined in example 3.3.1. The co- 
efficients for the required approximation are given, from equation (3.4.7), 
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Or 


Figure 3.4.3 The error associated with figure 3.4.2. 


as Co = 2.53213212 and c, = 1.130318, which gives the Chebyshev ap- 
proximation 


exp(x) & 1.266066 + 1.1303187; (x) = 1.266066 + 1.130318z 


which is close to the minimax approximation of degree 1. (Notice that 


the first term is $c.) 


Exercise 3.4.2 (Computational) Use the program APPROX with the 
option Method: Chebyshev to find approximations for the function 
exp(z), x € [-1, 1], of degree 1,2 and 3. Compare the computed approx- 
imations with the minimax approximations derived in the last section. 
(The resulting Chebyshev approximations are shown in figure 3.4.4.) 


Given the coefficients of a Chebyshev approximation we can com- 
pare it with the minimax polynomial. For example, the minimax cubic 
approximation for exp(z), z in [—1, 1], is given by 


p3(z) = 0.9945858 + 0.9956685z + 0.5429732z? + 0.17953272° 


which has minimax error €3 = 0.0055216. The corresponding Chebyshev 
approximation is 


P3(x) = 1.2660667p + 1.1303187; + 0.2714957> + 0.044337T3 
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Figure 3.4.4 Chebyshev approximations for exp(z). 


which gives the equivalent polynomial 
P3(x) = 0.994573 + 0.9973072 + 0.5429862? + 0.1773482°, 


and this has maximum error 0.006065. At this point we emphasize that 
there is no need to convert a Chebyshev series in this way; in fact to 
do so may lead to significant errors. The reason for this is the size of 
the coefficients in the corresponding series; the smaller the coefficients 
the better. Notice that for this example the coefficients of the higher 
powers of the Chebyshev series are much smaller than those of the min- 
imax approximation. It is better to evaluate a given Chebyshev series 
by making use of the recurrence relation (3.4.2) in order to minimize 
induced errors. In general given a Chebyshev series 


it can be evaluated as follows. For a given value of z define un(xz) = en 


Unie) eC ieats asl trys 
2 Crt t yd, aioe. Le (34.9) 


then 


\| 
oO 
pos 
— 

8 
— 

\| 

Sh. 
— 
8 
— 


tig( 2) = 5C0 + ru;(r) — u2(z) 
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Notice that (3.4.9) is a second-order recurrence relation, for any fixed 
value of x, and is relatively stable provided |z| < 1 which is precisely the 
range of values in which we are interested. (This technique for evaluating 
Chebyshev series can be compared with nested multiplication for power 
series.) In order to show this we rearrange the equations (3.4.9) to give 


Ch ae tna), 
Caer = tneila) — 2o ee 
cj; = uj(z)— 2ru;+1(2) + uj42(z), ) = 2 oe 
5C0 = uo(x) — ru;(z) + uda(z). 


Now recall that 


(Co) = Ben leer a $Co 
= un(2)(Tn(2) — 22Tn—1(2) + Tr-2(2)) + 
tinn1(2)(Ina1(@) — 2eI noe) in -al2)) ot 
u1(T, (2) — xTo(x)) + uo(z)To(z). 


Finally, we use the recurrence relation (3.4.2), in which case all the 
terms except the last on the right-hand side vanish; this last term is uo 
since Ty = 1. This technique would not be used to produce a explicit 
expression for a Chebyshev series as a power series but to evaluate the 
Chebyshev series at a particular point. 


Example 3.4.3 (Calculator) 
Evaluate the Chebyshev series 


1.2660667 9 + 1.1303187; + 0.2714957> + 0.044337T3 


at the point z = t We set u3 = 0.044337 and then produce 


ug = 0.271493+4+2 x 5 x u3 = 0.271493 + 0.044337 = 0.315832 

1 = ¢; 42x : X ug — ug = 1.130318 + 0.315832 — 0.044337 
seicd0T Sis 

Gm = 5C0 + sul — ug = 1.6511405. 


(Recall that co = 1.266066 from equation (3.4.7).) 


In general the determination of the Chebyshev coefficients is not as 
simple as demonstrated in example 3.4.1 and it may be necessary to 
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resort to numerical integration in order to evaluate them. (This is what 
the option Chebyshev in the program APPROX does.) However, the 
Chebyshev expansion (3.4.5) has one other bonus. Clearly, given any 
other set of polynomials we could expand any suitable function in terms 
of a similar expansion in which case we would still need to determine 
the relevant coefficients. Amongst all’ such expansions it can be shown 
that the coefficients in the Chebyshev series decrease more rapidly than 
any other expansion. Furthermore, given an expansion in terms of the 
first n such polynomials the error in truncating the series after n terms, 
1.e. 
CAT af Crit Ar oe 


is closely bounded above by the single term |cnT7,|, provided the coef- 
ficients tend to zero sufficiently rapidly. (See Fox and Parker (1972).) 
However, |T7,(x)| < 1 and so we can bound the truncation error by the 
single term |c,|. Therefore, if we require an approximation which has 
an overall error bound ¢€ we determine successive Chebyshev coefficients 
until we have |c,| < € and this will give the required polynomial. 


Example 3.4.4 (Computational) 
The cubic Chebyshev approximation for exp(x), z € [—1, 1] is given by 


P3(x) = 1.26606675 + 1.130318T, + 0.271493T2 + 0.04433773 


where the next coefficient is given by cq = 0.005474. The corresponding 
cubic minimax is given by 


p3(x) = 0.9945858 + 0.99566852 + 0.542973227 + 0.17953272°, 


which has minimax error €3 = 0.0055216. From this we can see that c4 
is a reasonable estimate for €3. 


In many cases Chebyshev approximations will be very close to the 
minimax polynomial. To illustrate this we list in table 3.4.1 various 
functions and compare the minimax error, corresponding Chebyshev 
maximum error and also the size of the first coefficient ignored in the 
Chebyshev approximation”. If the Chebyshev approximation is not suf- 
ficiently close to the minimax approximation then we simply take it as a 


1Strictly speaking amongst all expansions in terms of ultra-spherical polynomials 
for which the associated weight function is w(r) = (1—2x)~?, p= 4 gives Chebyshev 
polynomials, p = 0 gives Legendre polynomials. 

2The values given are those calculated by the software. If exact coefficients are 
supplied then different entries may be obtained. 
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Table 3.4.1 Minimax and Chebyshev maximum absolute errors. 


Taylor series Minimax Chebyshev Coefficient 


Order error error error neglected 

3 0.04744 0.005529 0.006063 0.005474 

exp(z) 4 0.00656 0.000547 0.000587 0.000522 
z€[-1, 1] 5 0.00292 0.000045 0.000064 0.000046 
2 0.28737 0.029830 0.033308 -0.029437 

log (1 +27) 4 0.15340 0.003424 0.003887 0.003198 
€[-1, 1] 6 0.07556 0.000444 0.000689 -0.000432 
sin(z) 3 0.00897 0.000499 0.000503 0.000521 
zé(-1, 1] 5 0.00009 0.000003 0.000003 -0.000003 


first approximation in the Remes algorithm. Alternatively, we can start 
the Remes algorithm by constructing a polynomial approximation using 
the points where the maximum errors occur at the extrema of the re- 
quired Chebyshev polynomial. For example, if we require a polynomial 
approximation of degree n, then we know that the error function must 
have n + 2 alternating minima/maxima; we take as our first approxi- 


ats) f= 0,1,...,.041, 
which are the n + 2 extrema of 7,4;(xz). This choice of points should 
produce an approximation which is close to the minimax. 


mation for the extrema the points 2; = cos 


Exercise 3.4.3 (Analytical) As an alternative to the Chebyshev polyno- 
mials, given on page 68, we can define the Modified Chebyshev polyno- 
mials T(x), for 0 < x < 1, according to 


CA Ort Rp Sai. 


Derive the first five Modified Chebyshev polynomials and sketch the 
corresponding curves. Show that 


nti(2) = 2(2@ — 1)T) (x) — Tr_, (2). 


The best HRER Dine Chebyshev approximation for r* +2? over the range 
[—1, 1] is x? + 22 which has a maximum error 0.25, i.e. 


eo 427 = 2? 4 $2 + 1T3(z). 
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Show that the best Chebyshev quadratic approximation for this function 
over the range [0,1] is given by 


pie (LUM STs eh 6Tt) Hato pe 


Therefore, the best Scans approximation is 35(1- 182+ 8027) which 
has a maximum error of + 35) One eighth that of the approximation on the 
interval [—1, 1]. 


Exercise 3.4.4 (Analytical). Runge’s problem revisited. In example 2.1.3 
and exercise 2.1.4 the construction of a polynomial approximation for 
exp(x) was considered. It was demonstrated that if the interpolation 
points are equally spaced then as the number of points increases the 
interpolating polynomial, p(x), oscillates wildly. This would seem to 
imply that 
2 
aes! len ( a) | = pes |exp(—z~) — pn(x)| + co as n > ov. 

However, recall that in general the error between an interpolating poly- 
nomial p,(z) and a function, f(z), for which it is an approximation 
satisfies 


Is(2) — pa(2)| < 2 Tle-- 


where M, = maxz,<c<z, |f)(zx)| and z;,i=0,...,n, are the interpo- 
lation points selected. Show that by selecting 


6 


n 


(x — 2;) = 27" Thai (2) (3.4.10) 


ll 
° 


i 
the error function e,(z) satisfies 


Mn+1 
6s. 
lenl@)IS 5 Ty 


Therefore, provided M,,4; is finite for all n, en(z) + 0 as n — oo. The 
choice (3.4.10) looks rather strange but simply requires that the zeros 


of both sides of this equation are identical, i.e. we need to select the 
interpolation points to be the zeros of 7,4; which are 


241 . 
ni = cos (AU) FXO) 10 Stern: (3.4.11) 
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> 3.4.1 Polynomial economization 


Let us suppose that we are given a polynomial of degree n and ask is it 
possible to construct a polynomial of degree < n which is minimax? This 
apparently difficult problem is in fact very easy to solve using Chebyshev 
polynomials. In order to see this let 


Qn(rz) = ao taie+...4+ Git) banca 


and then note that 7,(z) = 2"~!x" + P,_1(r), where P,_1 is a polyno- 
mial of degree n — 1. Therefore, 


2” = 20-)(T, (x) — Py-1(z)) 
which gives 


Qalz) =) apfaiete..+an-in" (and oo) in( Peal) 
= dofaie+ tae” tant ey Patt) hee cme 


Therefore, the (n — 1)th degree polynomial 
dot dieu. 4 a, ee an20-") P,,_ (x) 


differs from Q, by a term a,2\'-”)T,, (x) which has exactly n + 1 equal 
and opposite minima/maxima of size a,2\'—”) and is therefore the min- 
imax approximation of this degree. 


Example 3.4.5 (Analytical) 
Construct a quadratic minimax approximation for the polynomial 


Q3(z) = 1ltat+2?+ 2°. 
The Chebyshev polynomial 73(z) = 4x° — 3z, therefore, 
x? oe +(T3(x) + 32) 
and 
Q3=ltete2?+2° = l+a2+z2* +4(T3(x) + 32) 
= 144242? + 173(z). 


Therefore, the function 1 + ig + 2x? is the required quadratic minimax 
approximation for Q3(z) since it will have precisely four equal and op- 
posite minima/maxima of size i. No other quadratic polynomial can 
have a smaller error. 
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Exercise 3.4.5 (Analytical) Construct a linear minimax approximation 
for the quadratic function 


Q2x(z) =1+2+4+ $a’ 


which is the Taylor series expansion of exp(x) about x = 0. Determine 
the minimax approximation error involved. 


> 3.4.2 Chebyshev minimax approximation 


In section 3.3 the Remes algorithm was used to construct the minimax 
approximation for a given function. This algorithm is relatively easy 
to implement as was demonstrated in exercise 3.3.2 but as the order of 
the approximating polynomial increases convergence may be slow and 
it is possible that the process becomes ill conditioned. In the previous 
section it was shown how to construct a Chebyshev approximation which 
is almost as good as the minimax and so it would seem sensible to 
combine the Remes algorithm together with Chebyshev polynomials; we 
will briefly investigate this further in this section. 

Let f(z) be a continuous function for which p,(x) is the minimax 
polynomial of degree n. The error function e,(x) = f(z) — pn(z) has 
exactly n+2 alternating minima and maxima of equal magnitude. How- 
ever, aS p,(x) is a polynomial of degree n we can express it in the form 


Pn(v)+ = ScoTo +17; (x) + coTo(z) +... + enTn(z). 


This form of expressing the minimax polynomial has several advantages 
over an expansion in simple powers of z, not the least of which is that the 
coefficients will be numerically smaller which means that the evaluation 
of this expression will be less prone to rounding errors. Any other linear 
combination of Chebyshev polynomials will have a larger maximum error 
bound than that given by the minimax error €, = max|en(z)|. Let 
us assume that we have obtained an approximation for the minimax 
polynomial in terms of a Chebyshev series and that we have determined 
(n + 2) points at which the associated error function has alternating 
minima and maxima, zo < 2; <... < 2n41. Now we determine a better 
polynomial approximation, n(x), as 


Qn(z) = 5 bo + b, Ti (zx) +.s.+ Baliye) 
by solving the (n + 2) linear equations 
f(2:) — dn(2i) = (=1) t+ 8: amd Uae Ae al (3.4.12) 
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to determine the coefficients bo, 61, ..., bn. The similarity between the 
equations (3.4.12) and (3.3.3) is striking but the former equations have 
the advantage of being far better conditioned as n increases. In order 
to demonstrate this procedure consider example 3.4.6. 


Figure 3.4.5 The error exp(x) — Q3(z). See equation (3.4.12). 


Example 3.4.6 (Computational) 


Let us construct the cubic Chebyshev minimax approximation for the 
function exp(z),-1 < 2 < 1. We begin by specifying five points as 
extrema for the error function 


e3(x) = exp(x) — colo — C1 Ty — CoT> — €3T3. 


We use the extrema of the Chebyshev polynomial T4(z) which are given 
bya, = "cos (fir), 2 = 0,1,2,3,4. (See page 76.) This gives the points 
Do te pe —5 2.259 = 0t5— > 2 and x4 = 1. (Compare these 
values with those used in example 3.3.3). We now solve the five linear 
equations 


exp(zi) — coTo(2i) — c1Ti (a; ) — c2T2(2;) — Cs ls( 2.) (-1)'E, 
yn HO 823 A 
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These equations can be rewritten as 


Lives 1 —1 1 Co exp(—1) 
LGV 2eee Oni ts/25,—1 C1 exp(—5V2) 
1 0 —1 0 1 C30 j= exp(0) 
1. =5V/2 0 if? —1 C3 exp(4/2) 
1 iT 1 1 1 E exp(1) 


and produce the cubic Chebyshev polynomial 
Q3(z) = 1.266066 + 1.130321T, + 0.271540T> + 0.04488073 (3.4.13) 


which has the error function shown in figure 3.4.5. This function has 
maximum absolute value 0.00583 but the corresponding value of FE is 
0.005474. We see that this cubic is not minimax but we can repeat this 
process to improve it. We can also compare this cubic with the first four 
terms of the Chebyshev series for exp(z) given by 


exp(x) + 1.266066 + 1.1303187, + 0.2714957> + 0.04433773 


which has a maximum error of 0.00607. 


Exercise 3.4.6 (Computational) Write out the equations in the next step 
of the Remes algorithm for determining the cubic Chebyshev minimax 
approximation for exp(z). Solve the corresponding linear equations and 
investigate whether the resulting cubic is minimax. 


> 3.5 Orthogonal functions 


The Chebyshev polynomials introduced in section 3.4 are extremely use- 
ful when an approximation to a given function is required. This is be- 
cause they possess many useful properties, one of the most important 
of these being orthogonality. By this we mean that for a set of polyno- 
mial functions, ¢n(z), -1< c<1,n=0,1,..., there exists a weight 
function w(x) > 0 such that 


1 
{h w(t) On (2)bm(x) de = Mndmn, 


where M,, is a constant and bmn is 1 if n = m and 0 otherwise. We 
could take any other set of linearly independent polynomials, ¢,(2), 
x € [—1,1], and consider expansions similar to (3.4.1) or (3.4.5). (The 
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condition of linear independence is needed to ensure that no one member 
of the set can be expressed in terms of a linear combination of the 
rest and would, therefore, be redundant.) Ignoring the problems of 
convergence of infinite series, let us for the moment write 


f(z) = codo + 191 + cogo+.... 


However, if the functions ¢9, ¢1, ...are orthogonal we can determine 
the coefficients c, as 


Cn = M,, H fw x) f(z) onl x) dz 
Example 3.5.1 


The following are sets of orthogenal polynomials 
(i) Chebyshev polynomials. Ta(x), x €[—1, 1], w(x) = (1 — 2?)~/?. 


To(@ mel asia) =<; 


Dee) ee a 2 ee nd 
(ii) Legendre polynomials. P,(x), c € [—1,1], w(x) = 1. 
1 ea) eel eA Ee ee 
(n-- 1) P44(2) + abe) (2n 4) 2 (a): 
Further details of these polynomial functions and their properties can 
be found in Fox and Parker (1972) or Churchill (1963). 


The idea of orthogonal polynomials can be generalized further to sets 
of orthogonal functions. 


Example 3.5.2 

The following are sets of orthogonal functions 
(i) Exponential functions exp(inz), i = \/(—1), w(z) = 1. 
(ii) Fourier series. 


da(a) = { Sn) (a) = 1 


cos(n7z) 


Exercise 3.5.1 (Analytical) Show that the sets of functions given in ex- 
amples 3.5.1 and 3.5.2 are orthogonal. In each case, assuming that there 
exists an expansion of the form, 


n=oco 


f(z) = 9m cro, x), (3.5.1) 


n=0 
determine an expression for the coefficients c,, n= 0,1,... 


FOURIER SERIES EXPANSIONS 83 


Whenever we express a function in terms of an infinite series such 
as (3.5.1) we must always consider whether the series representation is 
valid, i.e. does the series converge and if so does it converge to the 
function we require. If the function is continuous then there are few 
problems, however, even when the function is discontinuous it is some- 
times possible to determine a suitable expansion. We will consider this 
problem in the next section. 


> 3.6 Fourier series expansions 


The idea of orthogonal functions was briefly introduced in the last sec- 
tion, we can demonstrate how useful this notion can be. Consider the 
following expansion in terms of the orthogonal functions cos(n7z), i.e. 


f(x) =cos(2mz) + }cos(6mxz) + 4 cos(10mz) +... 


+ iCos(2( 20 11) Aa)) toot (3.0.1) 


1 
(2n — 1) 


For the moment let us ignore the problem of whether this infinite series 
converges and try instead to plot it. Of course we cannot add in all the 
terms but for the moment let us see what happens if we include the first 
four, five or perhaps six terms. Now we have a finite expression which 
we can plot. The result is shown in figure 3.6.1. 


At first sight it appears that as we include more and more terms 
the sum of this series of continuous functions appears to be a function 
which is not differentiable at various points. Conversely, we seem to 
be able to write the sawtooth waveform as a sum of eminently smooth 
functions. We can now ask whether this is always the case. However, 
before we answer this question let us take a closer look at the graph 
drawn in figure 3.6.1, in particular near the peaks. By plotting in more 
detail we see that the graph in fact does not change direction sharply but 
does so in a smooth way (see figure 3.6.2). However as we include more 
and more terms the magnification required to see this smooth change 
becomes greater and greater and to all practical purposes the sawtooth 
wave is generated. 


Let us now turn to a general function f(z), |z| < 1, and examine 
the conditions under which it is possible to write this function as an 
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Figure 3.6.1 The cosine series (3.6.1). 


-0.95 
=11 (0/0) 
= 1,05 
= le ©) 
= Lie 
=. 20) 


SUAS) 


0.38 0.60 


Figure 3.6.2 Enlargement of figure 3.6.1 near z = OD =e ss 
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expansion of the form 


do n=0o 
Ec) aa 2 ( a, cos(ntz)+b,sin(nrz)), |x| <1, (3.6.2) 


where the ‘ denotes that the first coefficient is given by $40. The expan- 
sion (3.6.2) is called the Fourier series expansion of f(x). If we are given 
a function f(x) fora < x < b then we can always map the function to the 
interval |x| < 1 by a suitable linear transformation. Alternatively, if the 
function is defined on an infinite interval then since the sine and cosine 
functions are periodic, i.e. sin(z+27) = sin(z) and cos(zr+27) = cos(z), 
we shall require that the function f(z) is also periodic. Therefore, with- 
out loss of generality we shall only consider functions which are defined 
for |x| < 1. Next we consider the notion of sectional continuity. 


Definition 3.6.1. Let f(x) be a function which is continuous at all points 
of a finite interval a < x < b except possibly at a set of points a < 2, < 

.< 2, <b. If f has a finite limit as z approaches the ends of each 
interval (x;, 2441), 7 = 1,...,p—1, from its interior then f is said to be 
sectionally continuous. 


Definition 3.6.2. A function is said to have a right- (left-) hand derivative 


at a point if 
an £244) = F(@) 
im h 


h—0 


exists, where h tends to zero through negative (positive) values. 


Figure 3.6.3 Sectionally continuous functions. 
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Example 3.6.1 
The functions shown in figure 3.6.3 are simple examples of sectionally 
continuous functions. 


Theorem 3.6.1. Let f be asectionally continuous function on the interval 
(—1,1) then the Fourier series 


549 + (an cos(nmx) + by sin(nz)), 


il 


where 


Q 
3 
i} 


1 
i! f(z) cos(naz)dzjain = sie 2a 
mt 


1 
be ii flayon(ara da. er 
-1 
converges to the value 


aLf(z + 0) + f(z — 0)] 


at every point where f has a right- and left-hand derivative. (By f(z+0) 
we shall mean the limit of f at x through values which approach z from 
above or below.) 


Example 3.6.2 (Computational) 
Load the program APPROX and select the Method: Fourier. We will 
now investigate the Fourier series which correspond to the function 


Il, alice 7B Se 10), 
se) = { 0, 0 — eal. 


The function editor can be used to input this function piece by piece. To 
return to the function editor press ESCAPE and then when the prompt 
shows Action: Continue press RETURN which should produce the op- 
tion Function: as set. (Notice that the current function and range is 
shown in the bottom window.) By pressing any key except RETURN 
cycle through the options available at this time until it reads Function: 
inspect/edit and then press RETURN again. This will display another 
prompt which you will need to alter to Change f(x) before pressing RE- 
TURN. A flashing cursor will appear in the bottom window and we can 
type in the value 1 which now becomes the current function. The next 
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stage is to split the range using the prompt Split section, recall that we 
need to split the range at x = 0. The bottom window will now show 
that we have the function f(z) = 1 for negative values of z as the first of 
two sections. Next we shall amend the second section using the prompt 
Next section and then cycling until the prompt reads Change f(x), in 
which case we can set the function to f(z) = 0. Finally, we can plot the 
current function using the option Plot function. (There are several other 
options available at this level which can be used to move between the 
various sections of the function, edit or delete sections, or add other sec- 
tions if required.) At this point the required function should be plotted; 
notice in particular the form of the curve near the discontinuous point 
xz = 0. This is caused by using only a finite number of points when plot- 
ting. Now quit the function editor and return to the prompt Method: 
Fourier and then press RETURN again. Again at this point there are 
various options available which include supplying the coefficients manu- 
ally, reading them from file, storing the current coefficients in a file and 
calculating the coefficients corresponding to the current function. This 
last option requires the evaluation of equation (3.6.2) and in general this 
is carried out numerically. Therefore, accept the option calculate coeffs 
and then set the upper limit on the number of terms to four. After a 
short pause, whilst the coefficients are being calculated, the option will 
return to Option: calculate coeffs. Again cycle through the possibil- 
ities available at this point and plot the function, approximation and 
error. Increase the number of terms included in the series and compare 
the approximations which are determined. Notice that each approxima- 
tion takes the value 5 at. the origin which agrees with the value given 
by theorem 3.6.1. How do you account for the behaviour of the Fourier 


series at +1? 


Erercise 3.6.1 (Analytical) Determine analytically the Fourier coeffi- 
cients for the function given in example 3.6.2 and compare them with 
those determined by the program APPROX. 


Exercise 3.6.2 (Analytical) Show that the function f(z) = Jel, |2l <1 
does not have a one-sided derivative at the point rz = 0, but that this 
function does have a convergent Fourier series on the interval |r| < 
1 which is valid at zx = 0. This demonstrates that the existence of 
one-sided derivatives is not a necessary condition for the existence and 
convergence of Fourier series. 


Exercise 3.6.3 (Computational/Analytical) Use APPROX to determine 
the Fourier series expansion of exp(zx), |z| < 1. Compare the result with 
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corresponding polynomial approximations. Investigate the behaviour of 
the Fourier series which is obtained for values |x| > 1. 


Exercise 3.6.4 (Computational) Determine the Fourier series associated 
with the function |sin(7z)|, |x| < 1. 


Exercise 3.6.5 (Computational) Load the program APPROX and then 
input the file RAMP. Plot the function which is specified by this file 
and then determine its Fourier series. Compare the coefficients obtained 
with the exact analytical results. Compare the function and its Fourier 
approximation as the number of harmonics increases. 


> 3.7 Rational approximation 


In previous sections we have seen how it is possible to determine approx- 
imations for continuous functions in terms of polynomial expansions and 
sectionally continuous periodic functions in terms of Fourier series. How- 
ever, if we are given a function which has an apparent singularity neither 
of these methods can give a satisfactory approximation. Consider the 
following example. 


Example 3.7.1 (Computational) 

Load the program APPROX and determine a cubic approximation for 
the function f(x) = (2? — 3x + 2)7!, |z| < 1.5. Using the interpolation 
points —1.5, —1, 0, 1.5 produces the cubic shown in figure 3.7.1. The 
poor accuracy is due to the behaviour of f near zx = 1. This demonstrates 
that attempting to approximate a function which has a singularity by 
continuous functions cannot lead to a satisfactory result. Let us recon- 
sider the function f(z) = (2? —3z+2)7! and set F(x) = (x—1)f(z). If 
we now determine a polynomial approximation for F(z) then we obtain 
a very satisfactory minimax approximation. This suggests that the cor- 
rect way of approximating the function f(z) is to seek approximations 
of the form 

Pa(2) 

Qm(z)’ 


where P,(z) and Q(z) are polynomials of degree n and m respectively. 
Clearly, we could multiply P, and Qm by any factor and so by conven- 
tion we assume that the coefficient of r™ in Qn is 1. The resulting 
approximation is denoted by Ram(z). 


f(z) % 
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Figure 3.7.1 Cubic approximation for (x? — 3x + 2)7}, |r| < 1.5. 


Example 3.7.2 (Analytical) 
In order to illustrate rational approximation let us determine R, ;(x) for 
the function exp(x), |z| < 1. Let the rational approximation R,; be 


ao + ayz 


fondle bo +z 


then we need to determine values for ag, a; and bg. As with polynomial 
interpolation we can only make R, and exp(z) agree at a finite set of 
points and since we have three coefficients to select we ensure that 


CX a, ee era Oates clase 


For example, let us use the tabulation points z; = —1l, z2 = 0 and 
z3 = 1; then 


at x = —1 (bo — 1) exp(—1) = ap — ay, 
Aye = : bop = ao, 
Avera | : (bp + 1) exp(1) = ao + ay. 


Solving these equations produces a; = —1, ao = bp = —2.1639534. This 
gives the following rational approximation for exp(z): 

_ —2.1639534 — z 

~ —2.1639534 + 2 


which is shown in figure 3.7.2. Notice that this function is singular at 
xz = —2.1639534 but this is outside the region where the approximation 


iy ,1(2) 


90 THE APPROXIMATION OF FUNCTIONS 


is required. The error function corresponding to this approximation, 
e1(x) = exp(z) — R1,1(z), is shown in figure 3.7.3 and has an error 
bound of 0.0575. This error bound can be compared with the linear 
minimax error bound of 0.2788 and cubic minimax error bound of 0.0552. 


=i 0 1 


Figure 3.7.2 Rational approximation R;,1(x), for exp(z), |z| < 1. 


Errorx20 < 


Figure 3.7.3 The error Ri 1(x) — exp(z), |z| < 1. 


In general if we seek a rational approximation, Rym(z), for a given 
function f(z) then we can match the function and the approximation at 
n+m-+ 1 points since this is the number of coefficients to be determined, 
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1.e. we select points z;,1=0,...,(m +m), and then solve 
Or Zi) p23) ln (23) eh = 0,25-,(n +m); 


Since the coefficient of r™ in Q,, is unity we can rewrite these equations 


a; x! am bj a! f(z) Sac n a dl Gy ya al aed (pe S10 bea Bits) 


This gives a system of (n + m+ 1) linear equations in the unknown 
coefficients which can be solved by Guassian elimination with partial 
pivoting, if necessary. 


Exercise 3.7.1 (Computational) Use the program APPROX with the 
option Method: Rational to determine the R22(x) approximation for 
exp(z), |z| < 1, at equally spaced points. Compare the results obtained 
by using other points. Also compare the results with R3, and Ry 3 
approximations for this function. 


Exercise 3.7.2 (Computational) Use the program APPROX to con- 
struct rational approximations, Rmn(z), for exp(z), |z| < 2, for which 
m and n are less than 4. Compare them with the cubic minimax ap- 
proximation. 


As with polynomial approximation of functions based on interpola- 
tion at a set of points, rational approximation can suffer from induced 
ill conditioning. This can be overcome to some extent by using ratio- 
nal Chebyshev approximations and further details are given by Fox and 
Parker (1972). However, there are other problems which can occur. 


(i) The equations (3.7.1) may turn out to be singular. 


Example 3.7.3 (Analytical) 
Let us consider a rational approximation, f,:(z), for the cubic function 


f(z) = 2? — 3x? + 42, 


which is determined by interpolation at r = 0, z= 1 and x = 2. The 
equations (3.7.1) become 


a) 0 ao 0 
| ay = il hey? 
12 -4 bo 8 


and these equations are inconsistent. 
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(ii) In example 3.7.1 we saw that the rational approximation Rj; for 
exp(zx), |x| < 1, determined using interpolation at —1,0 and 1 is 


—2.1639534 — z 
—2.1639534 + x 
which has a singularity at r = 2.1639534. Fortunately, this singularity 


is outside the region where the approximation is required, but this may 
not always be the case. 


Example 3.7.4 (Analytical) 
Consider the function 


Ry (2) — 


f(x) =-34 202-827, 0<2<2. 


Using the tabulation points 0, 1 and 2 produces the rational approxima- 
tion 


which has a singularity at x = ; where the original function is continu- 
ous. 


(111) As P, and Qm are polynomials they may have a common factor 
which can be cancelled out. 


Example 3.7.5 (Analytical) 
Consider the cubic 


f(z) = 0.52? — 42? + 10.52 —6 


for which 
iF yf (2) S35/(3)-= 3) f(s 


The rational function R2 determined from equation (3.7.1) is 


zr? —2Qr r(x —2 
Baal) =e 


= 2, ife#2. 


In order that the resulting function is continuous at x = 2 we need to 
define R21(2) = 2 which is inconsistent with the data provided. 


We have seen that it is posible that rational approximation may 
introduce singularities, but this possibility can be advantageous as the 
following simple exercise demonstrates. 
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Exercise 3.7.3 (Computational) The function f(r) = exp(x)/(1—2 cos(x)) 
has a singularity at z = im. Construct a rational approximation for f(z) 
by considering the function (x — $7) f(z). 


> 3.7.1 Rational minimax approximation 


Clearly, given a function f(x) there will be many rational approxima- 
tions Ram. We can now ask is it possible to find one rational approx- 
imation amongst this class which minimises the maximum error. It is 
possible to obtain results analogous to the Chebyshev equi-oscillation 
theorem and also de la Vallee Poussin but we refer the interested reader 
to Ralston (1965). Similarly, there is a procedure which corresponds to 
the Remes algorithm which determines the rational minimax approxima- 
tion iteratively and is usually called Maehly’s algorithm. This requires 
the solution of non-linear equations at each step. Unfortunately, due 
to the non-linearity of Maehly’s algorithm there may not be a unique 
solution, but, it is possible to show that there is only one solution which 
does not introduce additional singularities into the required rational ap- 
proximation. The determination of this unique rational approximation 
is very difficult and beyond the scope of this text. 


Execise 3.7.4 (Computational) Use the program APPROX to contruct 
a rational approximation for exp(z), |z| < 1, by interpolation at five 
equally spaced points. Examine the error between the computed ap- 
proximation and the function and adjust the interpolation points to try 
to reduce the maximum error. Recompute a new rational approximation 
based on the new interpolation points and examine the error function 
which is produced. 


> 3.7.2 Padé approximation 


Just as it is possible to consider a truncated Taylor series of a continuous 
function as an approximation, the equivalent expression in terms of a 
rational expansion is called a Padé approximation. To illustrate this 
consider the following example. 


Example 3.7.6 


We shall construct a Padé approximation for exp(z) which is of the form 
Lig Wb 1.€. 
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where the coefficients ao, a, and bo are selected so that the function 
and the approximation have as many derivatives as possible equal when 
xz = 0. Firstly, we re-arrange this expression to give 


(bo + z)exp(z) = a9 + air (3:02) 
then when zx = 0, 1.e. 
bp exp(0) = bo = ao. (3.7.3) 
Next we differentiate equation (3.7.2) to obtain 
(bo + 2) exp(z) + exp(z) = a1 
and then set z = 0 to produce 
bb +1 =a). (3.7.4) 


As we have three parameters at our disposal we can repeat this procedure 
once more to give 


(bo + x) exp(x) + exp(r) + exp(z) = 0 
at x = 0 which gives 6) = —2. Substituting back into equation (3.7.4) 


gives a; = —1 and into equation (3.7.3) to get a9 = —2. Therefore, we 
have the Padé approximation 


exp(z) aes 
p(z) + pate 
This approximation can be compared with the rational approximation 
obtained by interpolation at x = —1, re = 0 and x = 1 in example 3.7.2 
which was 
—2.1639534 — x 
Ha (a) = 


~ —2.1639534 + 2° 


The errors associated with these two rational approximations is shown 
in figure 3.7.4. 


Exercise 3.7.5 (Analytical and Computational) Compare the rational ap- 
proximation Rj 3, R22 and R31 determined by interpolation at equally 
spaced points in (—1,1) with the corresponding Padé approximation. 
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0.1 


Rational error 


Figure 3.7.4 Padé and rational, R; 1, errors for e”. 


> 3.8 Applications and conclusions 


We now return to the problem posed at the beginning of this chap- 
ter; how does a computer evaluate sin(1.4) or exp(—3.456)? Different 
manufacturers have different methods but we shall illustrate the basic 
principles which are used. In general no approximation will be valid for 
all values of the argument and it is usual to define a principal range and 
then reduce the evaluation for a particular argument to this range. For 
example, we could use the periodicity of the trigonometric function to 
produce a principal range of (0,27), in fact we can use other smaller 
principal ranges as we now demonstrate. One other point to notice is 
the heavy emphasis placed on binary arithmetic which is the basic form 
of computer arithmetic. 


Example 3.8.1 (Analytical) 

As a first example we shall illustrate how the function sin(x), x in radi- 
ans, is evaluated. (The actual implementation will differ from machine 
to machine but the basic idea is common to all machines.) 


1. If |z| < 1078 then return sin(x) 2. 


2. Reduce the argument by periodicity so that |z| < 57 and set 


t= 2x/ 7, then 


hua Ge plete ix eth vpatea 
n=0 
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cos(z af Be Tao 


n=0 


where T** is the nth modified Chebyshev polynomial and the coef- 
ficients A, and B,, are given in table 3.8.1. 


Table 3.8.1 Modified Chebyshev coefficients for sin(x) and cos(z). 


n Ay B, 
0 1.276278962 0.472001216 
1 —0.285261569 —0.499403258 
2 0.009118016 0.027992080 
3 —0.000136587 —0.000596695 
4 0.000001185 0.000006704 
5  —0.000000007 —0.000000047 


Accuracy. For |x| < 5m the absolute error in the evaluation of sin(z) 


does not exceed 0.25 x 1078 and for cos(z) it does not exceed 0.42 x 107". 


Example 3.8.2 (Analytical) 

Evaluate sin(1.4). The required argument lies within (—}7, $7) there- 
fore we set t = (2 x 1.4)/7 = 0.891267681. Next we use ihe recurrence 
relations 


TES) came lisee B14 ral 
Try) = 2021? — 1) Tt) — sea) 
Then 
T(t Tj (t) = 0.588716160, 


T3 (t?) = —0.306826567, Tz(t?) = —0.949983676, 
T3(t?) = —0.811714915, Tz(t?) = —0.005755699. 


from which we obtain sin(1.4) + 0.985449729 which has absolute error 
0.78 x 10~°. Similarly, cos(1.4) = 0.169967143 which has absolute error 
(83 65xa1 Oanue 


Exercise 3.8.1 (Calculator) Evaluate sin(1.0) using the algorithm de- 
scribed in example 3.8.1. 
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Example 3.8.3 (Analytical) 
The evaluation of exp(z). 


1. If e < —180 return a value 0. 


2. Divide z by log,(2) and write 
x = (4a — b— sc) log. (2) — d, 


where a,b and ¢ are integers with 0 < 6 < 3,0 < c¢ < 15 and 
0<d< i$ log.(2), then 


exp(z) = 1672—°2-°/16 exp(—d). 
3. Compute exp(—d) using a minimax approximation polynomial ap- 
proximation of degree 6. 
4. Multiply exp(—d) by 2~°/!® then halve the result 6 times. 


5. Multiply by 16% by adding the hexadecimal exponent a to the 
exponent of the previous operation to give exp(z). 


Accuracy. If —1 < x < 1 then the maximum relative error in exp(z) 1s 
less than 0.21 x 107!°, otherwise it is less than 0.43 x 107}°. 


Example 3.8.4 (Analytical) 
The natural logarithm is evaluated using a rational approximation. 


elie —aelereturn log.(2) =): 


2. Determine an integer n such that 2"-! < x < 2” and then put 
Pe ES es pee peta i ( £V/2)/(r + 5/2) and find 


l+u 
O$e fom 


Vv = 
20790 — 21545.27u2 + 4223.9187u4 
~ ™\ 10395 — 14237.635u2 + 4778.8377u4 — 230.41913u® 
then 


log.(x) © (n — $) log, 24+ v. 


Accuracy. The maximum absolute error in log,(x) is less than 1.02 x 


Ljaee: 
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Exercise 3.8.2 (Calculator/Computational) Use the algorithms given in 
examples 3.8.3 and 3.8.4 to evaluate exp(3.456) and log,(3.456). 


Not all standard functions use polynomial approximations. For ex- 
ample, the evaluation of the square root function uses an iterative ap- 
proach. 


Example 3.8.5 (Analytical) 
The evaluation of \/z. 


1. If z = 0 then return 0. 


2. Let x = 167?-4 f, where p is an integer, g = 0 or 1 and 3 < of eay 
then 


Vz = 16°47%/f. 
3. Take as a first approximation for /f 


_, (0.14605 + 1.8342f) 
2168 4a ee 
hie (0.98026 + f) 


4. Apply the Newton—Raphson iteration 


1 ak 
Yn41 = 5) Yn + — 
Un 


twice. The final iteration is taken as yo = $(Mi — ea) + = in order 
to minimize computational errors. 


Accuracy. The maximum relative error is 0.48 x 107°. 


Example 3.8.6 (Calculator) 

Let us determine /3. We can write 3 = 162-! 0.1875, therefore p = 1, 
q = 1 and f = 0.1875. Using these values we obtain yo = 1.67829862, 
Y¥1 = 1.73291159, yo = 1.73205102 and y3 = 1.73205081. This last value 
is correct to eight decimal places. 


Ezercise 3.8.3 (Calculator) Use the algorithm outlined in example 3.8.5 
to determine V1.76. Square the result and hence estimate the error. 
Compare the error with the bound which is given. 


APPLICATIONS AND CONCLUSIONS oe 


In this chapter we have seen how it is possible to construct both 
polynomial and rational approximations for a given function. It is dif- 
ficult to say which type of approximation will in general be the best; 
this will depend on the accuracy required and the frequency with which 
the approximation is going to be used. It is arguable that in most cases 
a Chebyshev approximation will be sufficient for most purposes. The 
determination of a minimax approximation requires considerable effort, 
but once determined only the coefficients in the relevant polynomial need 
to be stored and evaluation is particularly easy. 


> Chapter 4 


> Least Squares Approximation 


Of course the first thing to do was to make a grand survey of the country 
she was going to travel through. “It’s something like learning geography”, 
thought Alice, as she stood on tiptoe in hopes of being able to see a little 
farther. 

Lewis Carroll. 


> 4.1 Introduction 


In Chapter 2 the problem of fitting a function through a set of data points 
was considered. It was shown that as the number of points increases 
then polynomial approximation becomes more oscillatory. In order to 
overcome this problem the concept of spline interpolation was introduced 
whereby lower order polynomials were employed. In this chapter we 
shall consider an alternative approach; to do this consider the following 
example. 


Example 4.1.1 (Analytical) 
Newton’s law of cooling states that the rate at which a body loses heat 
is proportional to the difference between the temperature of the body 
and the surrounding ambient temperature Tp. Let T be the temperature 
at time t; then we can write 

dT 

ne —K(T — Tp), 
where K is some constant, which can be integrated to give 


T(t) = To + Aexp(—Kt). (4.1.1) 


100 
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In an experiment the ambient temperature was found to be 10°C, and 
the data in table 4.1.1 were observed; how can we find the most suitable 
value for the constant K? The values of T — 10 are plotted against ¢ in 
figure 4.1.1, as are the curves Ty) + Aexp(—Kt) for an assortment of A 
and K. The problem of minimizing the error between the data points 
and any one member of this family of curves is quite difficult. 


Table 4.1.1 Data obtained for example 4.1.1. 
t(min) T°C  log.(T — 10) 


1 65 4.007333 
2 48 3.637586 
3 30 2.995732 


80 
60 
10+956°** 
40 / 104+956°°* 
/ vs -0.5x 
10+90e 
10+956°™* 
20 


1 2 3 
Figure 4.1.1 Values from table 4.1.1 and 10 + Aexp(—K‘t). 


As an alternative we can rewrite equation (4.1.1) as 
log.(T — To) = B—- Kt, B = log, A, 


and plot the new variable y = log.(T — To) against t. We then have a 
linear relationship between y and t, and the slope of this line is related 
to the value K we require. For a particular choice of B and K it is 
unlikely that the selected line will pass through all the data points and 
so we write the error at a point t; as 


ej; = y(t;) = fic (4.1.2) 
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OS 1 2 3 


Figure 4.1.2 Plot of t against log,(T — 10) for table 4.1.1. 


where y(t;) is the value determined from the particular choice of the 
parameters and y; is the observed value at t;. In order to obtain the 
‘best possible’ fit we could attempt to select the parameters B and K 
to minimize the sum of squares, S, of the errors. In this case 


S = (B—K —4.007333)? +(B-—2K — 3.637586)? +(B—3K —2.995732). 

(4.1.3) 
This function has a turning point where 0S/0B = 0S/0K = 0 which 
gives two linear equations for K and B. Solving these equations produces 
the least squares approximation for the given data as 


log.(T’ — Ty) = 4.558485 — 0.5058005t, (4.1.4) 


or equivalently 


T(t) = 10 + 95.43878 exp(—0.5058005t). (4.1.5) 


The corresponding sum of squares, S, given by equation (4.1.3) is then 
1.234 x 10-*. Notice that (4.1.4) is a linear least squares fit for the 
function log.(T —Tp) but that (4.1.5) is not necessarily a non-linear least 
squares fit for T. We shall discuss linear least squares in some detail; 
non-linear least squares is far more difficult. A detailed discussion is 
given by Stoer and Burlisch (1980). 
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> 4.2 Linear least squares 


In the last section we used a simple demonstration which fitted a straight 
line by minimizing the sum of squares of the associated errors. We will 
now extend this approach to fitting a linear function to a general set 
of values and then generalize the method to fit a linear combination of 
functions. 

Let us suppose that we have been given a set of tabular values 
(z;,¥;), 7 =1,...,m, and that we wish to find a linear function, 


y=cotcrz, (4.2.1) 


which minimizes some overall measure of the difference between the 
tabular values y; and the corresponding y(xz;) as provided by the ap- 
proximating function. There is no unique way of providing a measure of 
the discrepancy between y; and y(a;). For example, we could use the 
sum of the magnitudes of the errors to compare different values of co 
and ¢}, 1.e. 


j=m ,=110 
DU v1 — 8 leo ciz, — | (4.2.2) 
j=0 y=0 


It is possible to determine values of co and c; which minimize this func- 
tional but it requires linear programming techniques beyond the scope 
of this text. (If such an approximation is found it is called an i, ap- 
proximation.) However, in practice it is easier to minimize the sum of 
squares of the errors, i.e. 


j=m 
(Came ia S > (co +12; — yj)’. (4.2.3) 


j=0 


This is called the least squares fit or an 2 approximation. We could 
select values for cp and c, at random and compute the value of S(c¢o,c1) 
but a much simpler way is to use the calculus. 


Exercise 4.2.1 (Computational) Load the program LSQ and change the 
prompt to Select data. Now change the current option to Data: from 
keys and clear the existing data. Select the prompt Data input format: 
tabular and input the entries in the first and third columns of table 4.1.1. 
(Recall that log.(65 — 10) can be input as LN(55).) The data points can 
be plotted on suitable axes by using the option Data: plot&/or edit 
which should produce a screen display like figure 4.1.2. (Recall that 
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these data can be stored in a file using the option Data: store on file 
but you will need a disc of your own which is not write protected.) 
Exit from this set of options by using the prompt Data: quit and when 
the option shows Method: Linear press RETURN. Select the prompt 
Option: interactive fitting and a default line will appear on the screen 
together with the current sum of squares of deviations (SSD) as given 
by equation (4.2.3). The position and slope of the current line can be 
changed by using the cursor keys. 


1. t | shift the current line with respect to the axis y = 0. 
2. «— rotate the current line about the middle of the screen. 
3. <> increase or decrease the increment by which the line is changed. 


As the line is moved the sum of squares of deviations (SSD) is changed 
and displayed in the top window. It is also possible to examine the 
current values of the coefficients in equation (4.2.3), co and cj, by se- 
lecting the prompt Option: edit/inspect coeff. Try to adjust the line 
to minimize the sum of squares of deviations. 


Exercise 4.2.1 illustrates some of the problems in finding a least 
squares approximation by trial and error; therefore, we will now con- 
sider a more direct approach. The functional S(co, c,), given by equa- 
tion (4.2.3), has turning points whenever 


Ose Las 
re Ge (4.2.4) 


therefore, necessary conditions for a minimum are 


j=m 
a; (co 4+ ¢12;— y;) =. 0, 
j=0 
=m 
(comer a; 1; t= iO} 
7=0 


Finally, these equations can be rewritten in the form 


j=m j=m j=m 
2 
Co ) Ij+cy ‘ LjYj; 
j=0 j=0 =O 


(4.2.5) 
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j=m Se 
(m+1)co+ ec) Daas = Ure 
j=0 j=0 


For the example in section 4.1 these equations reduce to the linear equa- 


tions 
6 14 co \ _ f 20.269701 
Sat cl pane 0.640651 i) e 
The solution of these equations is co = 4.558485 and c; = —0.5058005 


as in section 4.1. These coefficients produce the fitted linear function 
y = 4.558485 — 0.5058005t. 


Exercise 4.2.2 (Computational) Load the program LSQ and select the 
prompts Method: Linear and Option: calculate coeffs. The program 
will now determine the least squares linear fit by using equation (4.2.6) 
and then return to this prompt. It is possible to display the fitted func- 
tion by selecting the prompt Option: plot and examine the computed 
coefficients using the prompt Option: edit/inspect coeffs. Compare 
the computed least squares fit with the results of exercise 4.2.1. 


Exercise 4.2.3 (Computational) Repeat exercise 4.2.2 but investigate the 
sensitivity of the computed line with respect to small preturbations of 
the data points. Is this problem well conditioned? 


We could generalize the above argument to consider functions of the 
form 
y(z) =co+ciz Oot te cn ee 
but these can in turn be considered as a special case of 


i=n 


y(z) = codo +eidit+...+ enon = YS cidi(z), (4.2.6) 
1=0 
where ¢;(x), i = 0,...,n, are any suitable basis functions. It is clear that 


n <m but if n = m the problem reduces to constructing an interpolating 
polynomial through the data points. As before we attempt to determine 


coefficients c;, 7 = 0,...,m, in order to minimize the functional 
j=m 
Dl corey: ance SOG; — y(z;))?. (4.2.7) 
j=0 


At a turning point we shall have 
Jone 


5g 0 FSO ey, (4.2.8) 
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which produces a system of n+ 1 linear equations in the unknown coef- 


ficients c;, i= 0,...,n. The equations (4.2.7) may be rewritten as 
j=m t=n 2 
SliCo eae ye € - S| cibi(z; ) (4.2.9) 
j=0 i=0 


Therefore, differentiating with respect to c, produces 


7=m te 
ys bx (2; ) (» = Yate) = 0, kee 0, sles Ue (4.2.10) ; 
=O) Z—=0 


which can in turn be re-arranged to give 


i=n [j=m FoR 
SPA you len ds(zs) bicep nde (Ta lUjagak uly eee ae 
7=O0\e7—0 j=0 


Equivalently, we may rewrite these equations as 
Ac=b, 


where the components of the matrix A are given by 


iI 
asi = | $e(2;)4s(2;) 


j=0 
and those of the vector b by bk = yet bx(z;)y;. Notice that A 
is independent of the values y;, j = 0,...,m. The system of linear 


equations (4.2.11) are called normal equations and can be solved by 
Gaussian elimination without the need to partially pivot since the matrix 
of coefficients is symmetric and positive definite. 


Exercise 4.2.4 (Computational) Load the program LSQ for which the 
default data points are given in table 2.1.1. When the prompt shows 
Select method press RETURN and then use the option Method: Linear 
to determine a linear least squares approximation for these data points 
using the option interactive fitting and then plot it. Determine the 
least squares fit by using the option calculate coeffs and compare the 
computed results. Now select the option Method: Polynomial and fit 
polynomials of degree 2, 3 and 4. (Why is it not possible to fit higher 
order polynomials to these data points?) Modify the data points using 
the option edit/inspect points and observe the way the fitted function 
changes. 
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Linear 
1S} 


Quadratic 


=O f Wes) 


Quartic 


Figure 4.2.1 Linear and polynomial fitting for table 2.1.1. 


Exercise 4.2.5 (Computational) Repeat exercise 4.2.4 but fitting 

(i) Chebyshev polynomials of degree 2, 3 and 4. 

(ii) Fourier series of order 1 and 2. 

(iii) The functions 1, exp(x), exp(2z), .... 

For (iii) what is the maximum number of such terms which can be used? 


Exercise 4.2.6 (Computational) Repeat exercise 4.2.4 but try specifying 
any suitable set of basis functions of your choice. (Recall that they must 
be linearly independent.) 


Exercise 4.2.7 (Analytical) Show that the functional, S, given in equa- 
tion (4.2.7) may be written in the matrix form 


(Costes Cricail INC i y)' (Ne —y), 


where N is a rectangular matrix for which N;; = ¢;(x;), y and c are 
vectors whose components are yj, j = 0,1,...,m andc¢,i=0,...,n, 
respectively. By differentiating this form directly show that the normal 
equations become 

N’Nc=N’y, 


and deduce that these equations are identical to the equations (4.2.11). 
Furthermore, show that the matrix N is symmetric and positive definite. 


The matrix A in the normal equations (4.2.11) will in general be full 
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but if the basis functions ¢,(z), i = 0,...,n, satisfy 


j=m 
br (2; )os(2;) = 0, ifr~s (4.2.12) 


7=0 


then all the off-diagonal entries of A vanish. If this is the case then the 
functions ¢;(r), i = 0,...,n, are said to be orthogonal with respect to 
the points x;, j = 0,...,m. Under this condition the normal equations 
reduce to a diagonal system and the coefficients of (4.2.6) are given by 


«= pnd (25) Yi 
eee 


Given any set of linearly independent functions it is possible to generate 
another set which satisfies (4.2.12) but we shall consider an alternative 
approach in a later section. 


Exercise 4.2.8 (Computational) Load the program LSQ and input the 
data file SPLINE. Compare linear, quadratic and cubic least squares 
fits for these data. Explain why the fitted function is so poor. Repeat 
this exercise for the datafile SIN. 


as) patentee (4.2.13) 


tae 


> 4.3 Continuous least squares approximation 


In the last section we discussed the problem of fitting a function of the 
form 


y(z) = codo(x) +... + cndn(z), (4.3.1) 


to a set of discrete points. This idea can be generalised to consider the 
approximation of a given function f(x), |z| < 1, by expressions of this 
form. In order to do this we consider the functional 


o (Cope ena I. (1) - ¥ «i()) dx (4.3.2) 


7=0 


as a measure of the error between the function and its approximation. 
It is also possible to consider a weighted measure of the discrepancy by 
defining 


1 tI) 2 
S( Cone Cn = / w(x) Le - ¥ a)) dz, (4.3.3) 


1=0 
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where w(z) is a suitable weight function. We have assumed, without loss 
of generality, that the function f(x) is defined on the interval [—1, 1]; if 
not then we can always transform it to be so. As before we differentiate 
S(co,.--,¢n) with respect to c, to produce the corresponding normal 
equations, i.e. 


y ( i w(x)bz(r)¢i(z) az) ee / w(x)d;(z)f(z) dz, k=0,...,n. 
i (4.3.4) 


This produces a system of linear equations in the form Ac = b, where 


— fp w(x), (x)o;(x) dx 


b= [ e@ee@re) dz. 


However, if the basis function ¢;(z), 1 = 0,...,n, are orthogonal with 
respect to the weight function w(z), i.e. 


ip w(x)o;(z)¢;(z)dzx =0, iF), 


then the normal equations are particularly simple and the coefficients 
CG, t=2,..4, 2,,are given. by 


— 4.3.5 
Jo, w(x)di(2)? de ae 


Example 4.3.1 (Analytical) 

The Chebyshev polynomials form an orthogonal basis. (See section 
3.4 and example 3.5.1.) Hence, we can construct the Chebyshev least 
squares approximation for a function f(x), z € [—1,1] as follows. We 
consider an approximation in the first (n + 1) Chebyshev polynomials 
given by 


i=n 
f(z) ® S- eT; (2). 
7=0 
Then from equation (4.3.5) the least squares coefficients are given by 
1a (aaa tc) 2p He)Fe) , 


— ——— dr, g=- ————— dz. 
Oe ds f= 2) tJ J — 2%) 
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Exercise 4.3.1 (Analytical) Show that for the coefficients given in exam- 
ple 4.3.1 the functional S(co,¢1,..-,¢n) given in equation (4.3.3), has 
minimum value 


1 2 
/ Se ee 


=i JY(1 — 2?) 


Exercise 4.3.2 (Analytical) Show that the Fourier series approximation 


k=n 
f(z) & Bo a, So(a cos kx + by sin kx) 
2 kar 


with the Fourier coefficients defined by 


is 1 us ; 
= a f(x) coskz dz, b= = f(z)sin kz dz, 


vie -1r -—7 


is a least squares approximation for f(r), with unit weight, in |z| < 7. 
(The question of the convergence of the Fourier series as n — o0 will 
not be discussed here. For a full discussion see Churchill (1963).) 


> 4.4 Discrete least squares approximation 


In the last section we saw how it was possible to construct a continuous 
least squares approximation for a given function, however, to do so it 
was necessary to evaluate the coefficients by performing numerous inte- 
grations. We will now return to consider how it is possible to construct 
a discrete variant which avoids this problem. 

In section 4.2 we saw that given a set of points (2; Oe 


m, 


we could construct a linear least squares fit, at G 1 Pal 2) DY solvint 
the normal equations given by 


Ss Ss x(x; )i(x;) C= v Pe(x;)y;, k= Oe. 10 (4.4.1) 


1=0 p= 0 3=0 


but if the basis functions are orthogonal with respect to the tabular 
points, i.e. 


jan 
> o(zs)oi(zj) =0, i Fk, (4.4.2) 
3=0 
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then the coefficients in the least squares approximation are given by the 
simple form 


—_ Lj=o bil, )y; 
1 Sea = ) 
a0 Pi(25)? 
Therefore, given any function f(x) we can construct the least squares 
approximation to it if (m+1) points, (z;, f(z;)), j = 0,...,m, are known 
such that equation (4.4.2) is satisfied. Conversely, if we have m+ 1 
points which satisfy this condition then the interpolating polynomial 
constructed in this way will be the least squares approximation. 


In order to illustrate this approach consider the Chebyshev approxi- 
mation given by 


= ee ne (4.4.3) 


ix=n 


fziis \> a:Tj(2), lero. (4.4.4) 


s—0 


where the ’ denotes the fact that the first coefficient is halved. Then if 


the tabulation points z;, 7 = 0,...,m, are selected to be the zeros of 
intl) 1.€. 
2j7+1 
j= COs (F243 =), a=) anit (4.4.5) 
it is possible to show! that 
j=m 0, eee ty 
Teas) i, (as) 4 ost, as 0, (4.4.6) 
j=0 Mm. or =s=0; 


where the denotes that both the first and last terms in the summa- 
tion, for any r or s, are to be halved. In this case the coefficients in 
equation (4.4.4) are given by 


a= — Y" F(es)T(25). (4.4.7) 


This is the discrete Chebyshev least squares approximation. If n = m, 
i.e. the number of tabulation points is equal to the number of coefficients, 
then the ' in equation (4.4.4) is replaced by " for an exact fit. Notice 
that the tabulation points z; are not equally spaced. 


1See Cheney (1966). 
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Example 4.4.1 (Analytical) ' 
Let us consider the function exp(z), z € [—1,1]. The Chebyshev series 
which represents this function is given by 
exp(z) = 1.2660667) + 1.1303187; + 0.2714957> + 0.044337T3 
+ 0.005474T, + 0.00054375 +.... 


The Chebyshev discrete least squares approximation of degree 2 which 
is constructed by using least squares at the zeros of T, is given by 


1.26606675 + 1.1303157, + 0.27145079. 


The maximum difference between these functions is 5.045 x 107? and 
the sum of squares of deviations is 3.835 x 1073. 


We can compare the coefficients of a Chebyshev series and those of 
the discrete Chebyshev least squares as follows. Recall that the coeff- 
cients of the Chebyshev series, c;, are given by 


a= =f f(cos@) cosi6 dé 
m Jo 


and if we apply the composite trapezium rule to this integral with m+ 1 
points we find 


7 wt 
ae ~ f (cos 0%) cos i6,, 6, = Kae 
uh m 
k=0 
5) k=m 
Pie F(2%)Ti(te) = ai, 
Wh: k=0 


from (4.4.7), i.e. the Chebyshev series of f(z) is the Fourier series of 
f(cos(z)). 


Exercise 4.4.1 (Computational) Load the program LSQ and change the 
prompt to Select function. When the current option shows Function: 
disabled change it to enabled and press RETURN. Now change the 
option to Function: inspect /edit and change the function to e”, |z| < 1. 
Quit this section of the program and change the prompt to Select data. 
Use the option Data: standard grid to select a uniform grid of four 
points and then quit. The program should now return to the prompt 
Method: Linear; change the option to Polynomial and fit a polynomial 
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of degree 2. Compare the result with a Chebyshev least squares fit of 
the same degree. What do you deduce? Increase the number of points 
and repeat the comparison. Return to the option Select data and use 
the zeros of T,(r) as data points, construct a least squares fit for e”, 
|x| < 1, based on these points and then compare it with your results for 
equally spaced points. 


Exercise 4.4.2 (Computational) Construct an interpolating polynomial 
approximation for e”, |x| < 1, which passes through six uniformly spaced 
points. (Use INTERP). Compare the result with least squares polyno- 
mial fits for polynomials of degree 3, 4 and 5. Repeat this exercise using 
the zeros of Tg(x). What do you conclude? 


Exercise 4.4.3 (Computational) Construct the minimax cubic polyno- 
mial approximation for e”, |x| < 1, using the Remes algorithm in AP- 
PROX. Compare it with least squares Chebyshev approximations of 
order 3,4 and 5. Estimate the relative efficiency of determining these 
approximations. 


> Chapter 5 


> Integration 


‘Quid pro Quadrature?’ 


Although many functions are known to be integrable, that is to say they 
have an indefinite integrand, there may be no expression for the inte- 
grand in terms of elementary functions. Perhaps the simplest example 
of this is the integral 


1 
| exp(—2”) dz. 
0 


In this chapter we will see how it is possible to determine an estimate 
for such integrals. To do this we first find an approximation for the 
integrand, integrate the approximation analytically in order to estimate 
the integral and then estimate the overall error. This process is called 
quadrature. 


> 5.1 Simple quadrature 


Let f(x) be a continuous function; then we may use its Taylor series 
expansion about x = a to approximate it, i.e. 


z—a)?* 


Fle) = f(a) +(2—a)f'(a)+ == g(a) 4. P=" pia)... 


(5.1.1) 
Therefore, 


[ f(x) dz 


f(a) [2 —a}’ + f"(a) [fese") ie 


a 
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=a) a) f(a) 4(b a)? fi(a)-Het. 


By adding in sufficient terms, and providing (b — a) is sufficiently small, 
we might hope to estimate the integral. If we were to replace the expan- 
sion (5.1.1) by a Taylor series with remainder we might even be able to 
estimate the error but as we have seen the Taylor series expansion is not 
the best approximation available. As an alternative, for any set of points 


23, j =0,...,m, we can find the Lagrangian polynomial approximation 
g=n 

Pr(z) = )_ L;(2)f(z;), (5.1.2) 
j=0 


fF) ‘fk 
Sch CA ACA a 1S CaN Saree “Me - Cae 


where 6 depends on z and lies between zo and z,. Therefore, 


[ seas 7 [ (ool2) + enla)yte (5.1.4) 
‘ 37 (x)f(zj;)de+ | en(z)dz (5.1.5) 
1% [- 


= Sot (5.1.6) 

where 
aj = i L,(x) dz, (5.1.7) 
E, = FT (a — 23) 0+ (6) de (5.1.8) 


The approach now is to select the points 2; in order to minimize the 
error E,,. Notice that the quadrature formula (5.1.6) involves (2n + 2) 
parameters; the interpolation points z; and the coefficients aj, j = 
0,...,n, but that these values are not independent. For any arbitrary 
choice of interpolation points, usually called the abscissa, the coefficients, 
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which are called the weights, are determined from equation (5.1.7). As 
a result. we should be able to make the quadrature formula exact for 
polynomials of degree n since n + 1 points determine a polynomial of 
degree n exactly. This is confirmed by the fact that e, is proportional 
to f("+!) which vanishes for all polynomials of degree less than n + 1. 
Quadrature formulae of this type are called Newton—Cotes methods and 
will be considered in section 5.2. Alternatvely, we can generalize this 
approach and begin with the formula 


b j=n 
/ f(z) dx = )_ a; f(2j)+ En (5.1.9) 
a j=0 


and then select values for all the z;’s and aj;’s to make this a method of 
order higher than n + 1. Such formulae are called Gaussian quadrature 
methods and they will be considered in section 5.3. 


Example 5.1.1 (Computational) 

The program QUAD is intended as a demonstration program for sim- 
ple quadrature formulae, whereas INTEGR is intended for practical 
numerical quadrature. In order to illustrate how the choice and number 
of interpolation points influences the accuracy of numerical quadrature 
load the program QUAD for which the default problem is 


1 
| x dz. 
0 


Accept the prompts Function: as set and Select data by pressing RE- 
TURN. The prompt will now change to Option: integrate. Notice that 
in the graphical window three equally spaced points are placed on the 
default integrand z°. This is confirmed in the bottom window where 
N = 3 is displayed. Furthermore, as three points are available it is 
possible to construct an interpolating polynomial of degree 2 which is 
confirmed by n = 2 in the top window. Now press RETURN. The 
bottom screen will display the estimate for the integral which has been 
determined as the Current Estimate together with an estimate of the 
quadrature error as determined by approximating equation (5.1.8). It is 
possible to examine the interpolating polynomial which has been used 
to approximate the integrand by changing the prompt to Option: plot, 
pressing RETURN and then accepting the option Plot: approximation. 
In addition, the difference between the integrand and the approximat- 
ing polynomial can be highlighted using the option Plot: shade and the 
error can be displayed using the option Plot: error. 
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Erercise 5.1.1 (Computational) Load the program QUAD with the de- 


fault problem 
1 
| x dz 
0 


and by changing the prompt Select quadrature to Select data investi- 
gate how the number of points in a uniform mesh influences the accuracy 
of the computed estimate for the integral. Why are six points sufficient 
to give an exact answer? 


Example 5.1.2 (Computational) 
Load the program QUAD to compute an approximation for 


iL exp(—z”) dz 


1 


using a uniform grid of five points’. (Hint: Use the option Select func- 
tion to change the integrand.) Return to the option Select data and 
then select Data: plot &/or edit and change the interpolation points 
to be at x = 0,40.5,41. Quit the data editor and select the prompt 
Option: integrate. An estimate for the value of the integral will now 
be determined based on these points together with an estimate of the 
quadrature error. Notice that the previous estimate remains visible in 
order that you can compare it. 


Exercise 5.1.2 (Computational) Repeat example 5.1.2 for other sets of 
five interpolation points and try to determine the best possible choice, 
i.e. minimize the computed error estimate. 


> 5.2 Newton—Cotes methods 


Let us now return to consider Newton—Cotes formulae which are based 
on replacing the integrand by an interpolating polynomial. We begin by 
considering the simplest possible case where we replace the integrand, 
f(x), by a linear approximation p;(xz), a< 2 <b, Le. 


b—z 


Peralta (Ole —s(0), (5.2.1) 


1Some computers treat the expression EXP(—X?) as EXP((—X)*) whereas we 
require EXP(—(X)?). 
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which is the Lagrange polynomial through (a, f(a)) and (6, f(6)). There- 


fore, 


[toe = [tears FE 


fa) f° Okey prseniesys 
a fom a)det oy | \dz+ EF, 


= }(b—a)[f(a) + f(d)) + Fi. 


This is the simple trapezium rule. We may use equation (5.1.3) to 
determine the error, £, in this approximation as 


b b FN (Ox " b 
[ awa | FU e—ay(2-b) dz = OM) f (e-ay(e-0) de 


(5.2.2) 
where the mean value theorem for integration has been used to simplify 
the integral. Integrating the last term gives 


aoe 5 ) (0), (5.2.3) 


where a < 6 < b. The accuracy of the trapezium rule is proportional 
to the cube of the interval width, provided the integrand has a bounded 
second derivative, and if we require higher accuracy then it is necessary 
to include further interpolation points. 


Erercise 5.2.1 (Analytical) Show that by sub-dividing the interval (a, 6) 
into two sub-intervals such that z9 = a, 2; = $(a +b), 2 = 6b and then 
replacing the integrand by a Lagrangian polynomial approximation of 
degree 2 we obtain 


b T2 
: ft dz i! po(x)dx + BE» 
= ¢(b-—a) [f(zo) + 4f(a1) + f(z2)] + Eo, 


which is Simpson’s rule. Show that 


eo 
Ey = —— f’” Kt) K th, 
9 99 fO9): a<@<b 
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The trapezium rule will integrate any polynomial of degree < 1 exactly 
since the error is proportional to f’, whereas Simpson’s rule will inte- 
grate any cubic without error. In general a quadrature formula is said 
to have degree of precision p if it will integrate all polynomials of degree 
< p without error. The trapezium rule has degree of precision 1, Simp- 
son’s rule has degree of precision 3. By introducing further points we 
can increase the degree of precision further. 

Example 5.2.1 

The Newton—Cotes closed formulae for the integral ayy f(x) dz are given 
by the following. Divide the interval (a, b) into n sub-intervals such that 
Co 10 nls 6 + fit) ng Te) = (0 —a)/n then ie loys 
is given in table 5.2.1. Notice that for n even the degree of precision is 
n+ 1, but for n odd it is n and this is generally the case (see Isaacson 
and Keller (1966)). Accordingly, even formulae are usually preferred. 


Table 5.2.1 Closed Newton—Cotes formulae. 


n 


a Se SE ES eee 
=1 zh[f(zo) + f(z1)] s 
2 point The trapezium rule: p=1: E, = —h'f"(0) 


n=2 sl f(zo) + 4f(21) + f(z2)] 
3 point Simpson’s ard rule: p= 35 = —gah® f'"(8) 


n=3 h[f(zo) +3f(21) + 3f(22) + f(za)] ; 
4point Simpson’s 3th rule: p=3: E3= —Sh° f'"(8) 


Tea Ah[7f (xo) + 32f(21) + 12f (x2) + 32f(z3) + 7F(r4)) 
5 point Milne’s rule: p=5: Eq = —gSzh" f""(0) 


The formulae given in example 5.2.1 are termed closed since they 
use the end points of the interval of integration; however, it is useful to 
define formulae which do not have this property. 

Example 5.2.2 a 

The Newton—Cotes open formulae are given as follows. We divide the 
interval (a,b) into n + 2 intervals so that a = z_; and b = 2n}41, i.e. 
2j =at+(j+I)h, j = 0,...,n, with h = (b—a)/(n + 2); then open 
quadrature methods for ifs f(x) dz are given in table 5.2.2. 
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Table 5.2.2 Open Newton—Cotes formulae. 


n= Us e2n f(a) 
1 point The midpoint rule: p=1: Eo = 3h? f’"(8) 


n=1 Sh[f(z0) + f(21)] 
2 point pal beh 7 (8) 


v= 2 4h[2f(z0) — f(xi) + 2f(z2)] 
3 point p=3: E, = Hh fi*(6) 


Example 5.2.3 

Consider the integral ie he dz. It is not possible to apply the closed 
Newton-—Cotes formulae to this integral even though the integral is well 
defined because the integrand is undefined at z = 0. However, there is 
no such difficulty in applying open formulae which do not require the 
integrand at x = 0. 


Ve 


As the number of points increases we obtain formulae with higher 
precision but unfortunately the resulting methods become unsuitable. 
For example, the nine-point closed Newton—Cotes formula is given by 


/ i f(z) dx = =A4-[989f (x0) + 5888 f (x1) — 928f (x2) + 10496 f(z3) 
—4540f (24) + 10496 f(s) — 928f(z6) + 5888f (x7) + 989f(zs)]. 


This formula has large coefficients which alternate in sign and this can 
lead to large cancellation errors unless special care is taken. 


> 5.3. Gauss quadrature formulae 


We shall now attempt to find quadrature formulae of the form 


1 gets 
[ tee = olen ea, (5.3.1) 
a am 


which have degree of precision > n+ 1. Such formulae are not restricted 
to the interval (—1,1), since given any other interval we can simply map 
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it to this range. (The reason for using this particular range will become 
clear later.) We will now select the (2n+ 2) parameters, i.e. the weights, 
aj, and the abscissae z;, j = 0,...,n, to make the formula as accurate 
as possible. 


Example 5.3.1 (Analytical) 
For n = 1 equation (5.3.1) becomes 


1 
ifs f(z) dz = aof (zo) + aif (21) + Ej. (5.3.2) 


As we are looking for a method which has degree of precision > 2 the 
error E; must vanish if f(z) = 1. Therefore, 


[ bdr 2 =ao + 01. (5.3.3) 

Similarly, if f(x) = x then EF, = 0 provided 
[ ze 0=ao2%0 + a2}. (5.3.4) 

We can repeat this for f(z) = x? to obtain 
Pages haat (5.3.5) 


and for f(z) = x3 to obtain 
0 = apx3 +a42°. (5.3.6) 


If all the equations (5.3.3)—(5.3.6) are satisfied then the resulting quadra- 
ture formula will be exact for all possible cubics, i.e. it will have degree 
of precision 3. Notice that the solution of these equations would be iden- 
tical if we were to integrate from 1 to —1 and so ag = a and zp = —2} 
and these conditions enable us to determine 


1 
(of 6 hat A eda tan / 5, 


This gives the Gauss quadrature formula 


fsmwes(-f)1(h) + 


where E; is proportional to f'’(8), —1 < @ < 1, since this method has 
degree of precision 3. 
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Example 5.3.1 demonstrates the method of undetermined coefficients 
and clearly this approach can be generalized to the Gaussian quadrature 
formula given by equation (5.3.1). In equation (5.3.1) there are (2n + 2) 
parameters, (n + 1) abscissae z;, j = 0,...,n, and (n + 1) weights a;, 
j = 0,...,n. Therefore, we should be able to select the abscissae and 
weights to give a degree of precision of at least 2n+1. It is also clear that 
to use this method for the derivation of Gaussian quadrature formulae 
will lead to lengthy algebraic manipulation if n > 1. Fortunately, there 
is an alternative approach. 

Recall that a set of polynomials, ¢;(z), 7 = 0,.... -l1 < z < 1, are 
said to be orthogonal with respect to a weight function w(z) if 


1 
ik w(x); (x)ox(z) dx = 0, Veale 


We will now show that if we use the (n + 1) zeros of ¢n41(x) as the 
abscissa in equation (5.3.1) then we can obtain degree of precision 2n+1. 
To do this it is necessary to show that ¢,(z) has exactly n real zeros in 
the interval (—1,1). Firstly, ¢o is a constant C # 0, therefore, 


1 1 
C= ie w(r)bn(x)go(x) dx = C ib w(z)on(x) dx 


and since w(x) > 0 this implies that ¢,(r) has at least one change of 
sign in (—1,1) and hence has at least one zero. Recall that ¢, is a 
polynomial of degree n and by the Fundamental Theorem of Algebra 
has at most n real zeros. Let us assume that the number of real zeros 
of ¢, in (—1,1) is j < n, and that these zeros are at z;,7 = 1,...,], 
where —1 < 2, < 22 < ... < z < 1. With no loss in generality we 
may assume that én > 0 in (—1, 21), dn < 0 in (2, 22), and so on. Now 
define P(x) = (—1) J];=)(2 — z;); then P(x)¢,(x) > 0, and so 


/ w( 2) P (eon 0: (5.3.7) 


1 


However, P is a polynomial of degree j < n and so it can be expressed 
as a linear combination 


Aes) 


P(z)= cide, 


4=0 
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therefore, 


1 t=j 1 
i w(x)P(z)¢, dz = DEG; / w(x)¢o;¢, dz = 0 
a 2=0 ae 


which contradicts (5.3.7). Hence all n zeros of ¢,(z) lie within (—1, 1). 
We have already seen examples of orthogonal polynomials in Chapter 


3. In particular the Legendre polynomials are given by Po(x) = 1, 
Pie p= and 


nee DP) Pra = (2n 2 Pa (2) mn Pye) pH 1 <r 1 6 (5.3.8) 


which have weight w(z) = 1. From above we see that P,41(r) hasn+1 
real zeros in the interval (—1, 1), and we will now show that if we select 
these zeros as the abscissa we obtain a quadrature formula of the form 
(5.3.1) which has degree of precision 2n + 1. To show this let P(x) be 
any polynomial of degree < 2n +1, then there are polynomials Q(x) and 
R(x) of degree < n such that 


P(z) = Q(2)Pa4i(z) + R(2). 


Furthermore, we can express (Q(x) as a linear combination of Legendre 
polynomials of degree < n, 1.e. 


r) = SE d; (a) 
7=0 


hence 


[ie Q(z)Pa4i(a) de = ya [6 )Pn4i(x) dz = 0 


and so 


IE P(2) a= fl Q(z Pate) dot Rade =f R(x) dz. (5.3.9) 


Given any set of abscissae we can select weights a; such that 


i Fila) date 5 aj R(2;) (5.3.10) 
eon f=0 
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where the integration is exact since R(x) is a polynomial of degree less 
than n. Finally, notice that if the interpolation points are selected to be 
the zeros of P,4i(z) then 


P(2;) = Q(2;)Pn4i(2j) + R(z;) = 0+ R(xj) = R(z;) 


and so 


1 j=n 
ii PG\do ==) o;P@,) 
F = 


1 


with no error, i.e. the quadrature method has degree of precision 2n + 1. 
It is possible to show? that the error of such a formula is given by 


h2Pa2((n ae Hj ))2 


m5 (aE) (O(a) sah et a ma 


where h is the interval width. For the interval (—1,1) we have h = 2. 


Example 5.3.2 (Computational) 
In order to illustrate the effectiveness of using the zeros of the Legendre 
polynomials load the program QUAD and then estimate the integral 


1 
i xr°dz 
0 


using three uniformly spaced points. (This is the default problem for 
the program QUAD.) Press ESCAPE and change the prompt to Select 
data and then press return. Now use a standard grid of three points 
which are the zeros of the required Legendre polynomial. This produces 
an estimate for the integral of 0.1666669 compared with the exact value 
of é Recall that in example 5.1.1 it takes six equally spaced points to 
evaluate this integral exactly; here we have used only three. Notice also 
that the function which is used to approximate the integrand differs from 
2° by as much as 0.207 and yet the estimate for the integral is very good. 
The reason for this is given by equation (5.3.11). Using three points we 
have n = 2 and so the error term will involve f?'(0) which is identically 
zero for this integrand. (Notice that the program automatically deals 
with the problem of mapping the zeros of Legendre polynomials, usually 
defined on the interval +1, to a suitable range.) 


2 See Froberg (1987). 


GAUSS QUADRATURE FORMULAE 125 


Example 5.3.2 (Analytical) 
Let us find the Gauss—Legendre quadrature method based on two points, 
l.e. n = 1. The second Legendre polynomial is 


P(e) = a (occ — 1) 


which has zeros at zo = -/3 andsa = a3 these points are used as 
the abscissa. The weight ao is given by 


1 
Ch = Cah 
/ Beate 
-1 %0 — 21 


dz 


ao 


@-yp 
1 


Similarly a; = 1. This gives the quadrature formula 


1 
1 1 
f(a)dz x1 (-4z) +1() 
I, BB) 
which was obtained in example 5.3.1. Finally, the error is given from 
equation (5.3.11) as 4: f*"(0), |@| < 1. We see that this method which 
is obtained by approximating the integrand at two points has degree of 


precision 3. Recall that using two equally spaced points the maximum 
degree of precision is 1. 


Example 5.3.3 
Values for the weight and absicca of Gauss—-Legendre quadrature formu- 
lae are given in table 5.3.1. 


Example 5.3.4 (Calculator/Computational) 
Let us estimate the integral 


1 
| exp(—2z*) dx 
0 
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Table 5.3.1 Gauss—Legendre quadrature weights and abscissa. 


n Dy a; 
1 +0.5773502692 1.0000000000 
2 0 0.8888888889 
+0.7745966692 0.5555555556 
3 +0.3399810436 0.6521451549 
+0.8611363116 0.3478548451 
4 0 0.5688888889 
+0.5384693101 0.4786286705 
+0.9061798459 0.2369268851 


using Gauss—Legendre quadrature. The range of integration is (0, 1) 
which can be transformed to (—1,1) using the substitution t = 2z — 1 
but notice that 


1 1 
Ji exp(—27) dz = 2 f exp(—zx”) dz. 
0 -1 


Applying the Gaussian quadrature formula with n = 1, 1.e. two points, 
gives the result 0.7165315 whilst the exact value is 0.74682413.... This 
is rather disappointing but the reason for the low accuracy is clear if we 
sketch the integrand. From figure 5.3.1 the integrand has a relatively 
large peak at x = 0 which is ignored by Gaussian formulae which use an 
even number of points, i.e. n is odd. However, if we use the Gaussian for- 
mula with n = 4, 1.e. five points we obtain the estimate 0.746832, which 
is correct to four decimal places. Using the program QUAD applied to 
this integral directly, i.e. without changing the range of integration, and 
using the zeros of the second Legendre polynomial as abscissa produces 
an estimate for the integral of 0.7465947. Using the zeros of P4 produces 
an estimate of 0.7468245 whereas using four equally spaced points gives 
0.7469923. 


Exercise 5.3.1 (Computational) Load the program QUAD and estimate 
the integral 


1 
/ exp(—z”) dx 
0 
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-2 = 1 2 
Figure 5.3.1 The integrand exp(—z”), -l< 2 <1. 


using two, three, and four equally spaced points. In each case exam- 
ine the interpolating polynomial which is constructed and compare it 
with the integrand. Notice that the program computes an error bound 
between the integrand and the interpolating polynomial and also an es- 
timate of the quadrature formula. Repeat this exercise using the zeros 
of the corresponding Legendre polynomials. 


The Gauss—Legendre formulae were deduced from orthogonal poly- 
nomials which are defined only on the interval (—1, 1) and integrals over 
other ranges must first be reduced to this interval. However, it is fre- 
quently the case that other sets of orthogonal polynomials can be used 
to evaluate integrals defined on other intervals. 


Example 5.3.5 
Alternative weight functions. 
(i) Laguerre polynomials. (See table 5.3.2.) (eaw x) f(x) dz 


ih |e tOmane Oust COL, (2) =e” Sasa ia 


(ii) Hermite polynomials. (See table 5.3.3.) f°. w(x) f(x) dx 


w(x) = exp(—z”): -oo <x < ow. 


Hp(x) = (—1)" exp(2?)(d"/da")(exp(—2”)). 
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Table 5.3.2 Gauss—Laguerre weights and abscissa. 


n 


1 


Pe 
0.585786 
3.414213 


0.415775 
2.294280 
6.289945 


0.322547 
1.745746 
4.536620 
9.395071 


0.263560 
1.413403 
3.596426 
7.085810 
12.640801 


ie 
0.853553 
0.146447 


0.711093 
0.278518 
0.010389 


0.603154 
0.357419 
0.038888 
0.000539 


0.521756 
0.398667 
0.075942 
0.003612 
0.000023 


Table 5.3.3 Gauss—Hermite weights and abscissa. 


n 


1 


2 


t: 
+0.707107 


0 
+1.224745 


+0.524648 
+1.650680 


0 
+0.958572 
+2.0201829 


Qa; 


0.886227 


1.181636 
0.295409 


0.804914 
0.081313 


0.945309 
0.393619 
0.019953 
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> 5.4 Composite quadrature 


We have seen that the error associated with simple quadrature is propor- 
tional to a power of the interval width. It would seem sensible, therefore, 
to make the range of integration as small as possible. To do this we write 


b at+h a+2h b 
fi f(a)dz = f f(a)de + | f(r) dr+...4+ Keys 
a a at+h b-h 
(5.4.1) 
where h = (b—a)/N. We then replace the contribution made in each 
sub-interval by the corresponding quadrature formula. The result is 
called a composite quadrature formula. 


Example 5.4.1 (Composite trapezium rule) 
The simple trapezium rule applied to an interval (2;,2;4,) is given by 


eg f(x) dx = Sh[f (xi) + f(zit1)] + Ai, 


where h = (x; — 2441) and Ej} = —h?f"(0;), 2; < 0; < 2i41. If we set 
z; =a+th, with h = (b—a)/N, then equation (5.4.1) becomes 


b 
[tear = GMs) + fa +ml+ Bf 


PG +h) + f(a+2h)] + Ei + 


+5 hLf(b— h) + FO) + EN 


h[f(a)+2f(ath)+2f(at+2h)+...+2f(b—h)+f(6)] 
+Ecr, 


a 
2 


where the composite trapezium error, Ect, is given by 
t1=N-1 
he on NhP ey 
= — j= —-— 6). 
Wonte-ae de) agai foal Oe) maa FAB) 


*1=0 


(The mean value theorem for integrals has been used to replace 0;, 2 = 
0,...,N —1, by 0 with a < 0 < b. See Appendix A.) However, Nh = 
(b — a) and so 

Eor = —ph?(b- a) f"(8). (5.4.2) 
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Notice that the error of the composite trapezium rule is O(h?) compared 
with O(h®) of the simple trapezium rule, but that the h in the compos- 
ite rule can be made small by introducing more sub-intervals whereas 
with the latter it is the full interval of integration. By systematically 
increasing N we should be able to reduce the quadrature error. 


Example 5.4.2 (Computational) 
Consider the simple trapezium rule for the following integral: 


s[f(0) + f(1)) + Er 


1.8591409 + FE. 


1 
/ exp(x”) dz 
0 


We can bound the error, £;, involved in this procedure by using equa- 
tion (5.4.2), i.e. 


|E| = |-4(0)| < & 6exp(1) = $e? = 1.35914 


since |f’"(x)| = |(4x? + 2) exp(x?)| < 6. Therefore, 
1 
| exp(x”) dx = 1.8591409 + 1.35914. 
0 


We can improve upon this estimate for (is exp(z?) dz by using the com- 
posite trapezium rule. For example, loading the program INTEGR, 
for which the default problem is this inegral, we obtain the results given 
in table 5.4.1. Notice that the program gives an estimated error bound 
which is determined by estimating an upper bound for |f”| using fi- 
nite differences®. In all cases this appears to be an overestimate of the 
quadrature error which is also given in table 5.4.1. The actual error in 
determining this integral is given in the final column of this table. 


Notice that the errors in table 5.4.1 decrease by a factor of four when- 
ever the step size is halved, thus confirming that the error in the com- 
posite trapezium rule is proportional to h?. Clearly, as h tends to zero 
the value produced by the quadrature formula appears to converge and 
by systematically reducing h we should obtain an accurate result. Al- 
ternatively, we can estimate the number of intervals required to achieve 


3See Appendix B. 
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Table 5.4.1 Composite trapezium rule for ip exp(z”) dz. 


N Estimate Estimated error bound’ Error bound Error 


2 1.5715832 0.0958526 0.3397852 0.1089310 
4 1.4906789 0.0410165 0.0894463 0.0280267 
8 1.4697123 0.0143884 0.0212366 0.0070601 
16 1.4644203 0.0043409 0.0053091 0.0017681 


a given accuracy. For example, let us suppose that we require a final 
estimate which is correct to six decimal places, i.e. |Ect| < 0.5 x 107°. 
To do this we select h such that 


a ll phe " -6 
|Ecr| Ss ale cael lf | Sopa) 


or 


he 1 10-8 
€ 


h 0.6065 x 107°. 


IA 


From this we conclude that N = 1/h > 1649 intervals. Using this value 
for N gives 


1 
/ exp(x”) dx = 1.462651912 + 0.00000017. 
0 


By increasing N further we might expect to be able to increase the 
accuracy further but we are limited by the accuracy of the machine 
being used. To see this recall that any estimate for an integral which 
is computed by quadrature formulae of this form is determined by the 
addition of N terms and if each term is evaluated to a precision ¢ then 
the total has a maximum possible error Ne; as N becomes larger the 
accumulated error will dominate the quadrature error. 


Erercise 5.4.1 (Composite Simpson’s rule). Show that by replacing an 
integral by 


b 


[tears [ serars [ seyars...+ | f(x) dz 


+2h b-—2h 
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where h = (b— a)/2N, and then applying Simpson’s rule we obtain 


[ t@w= th[f(a)+4f(a+h) + 2f(a+ 2h) +4f(a + 3h)+ 
: ... + 2f(b— 2h) + 4f(b—h) + f(b)] + Ecs. (5.4.3) 


The error of this composite Simpson’s rule is Ecs = — rag ht*(b- a) f*” (0), 
Go! 


Erercise 5.4.2 (Computational) Load the program INTEGR and esti- 
mate the integral ie exp(—z”) dx using the option Trapezium. Use the 
error term Ecr to estimate the number of intervals required to deter- 
mine this integral correct to five decimal places. Repeat this exercise 
with the composite Simpson rule using the prompt Simpson's 1/3rd. 


It is clear that we can develop composite integration formulae for any 
of the higher order Newton—Cotes or Gauss methods. Notice that for 
hand calculation when using the Newton—Cotes formulae and doubling 
the number of points from N to 2N we can re-use previous calculations. 
For example with the composite trapezium rule with N intervals we have 


b 
| f(x) dz » T(h) = 1ALf(a) + 2f(a +h) +... 42f(b—h) + FO) 
and for 2N intervals 


T(5h) = $h[f(a)+2f (a+ 4th) +2f(a+h) 
+2f (a+ $h) +... 
+2f(b—h)+2f (b— 4h) + f(0)] 
= 17(h)+ 
sh[f (a+ 5h) +f (a+ 3h) +...4+f (b-4h)]. 


Therefore, to find T($h) we need only N more function evaluations 
than when finding T(h). This is generally the case with Newton—Cotes 
formulae. The Gaussian quadrature formula with 2N points has no 
abscissa in common with the N-point formula. The composite trapezium 
method has another advantage; the weights at each abscissa remain the 
same as the number of intervals is doubled but the same is not true with 
higher order methods. 
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Example 5.4.3 (Analytical) 
Consider Simpson’s rule applied with two intervals to produce 


b 
ij f(z) dz ~ Eh[f(a) + 4f(a th) + F(d), 


where h = $(b— a). If we double the number of sub-intervals then 


[ f(a)dzm ch [fa)+4y (at5h) +2setn+4s (at5h) +50) 


and we see that the coefficient of f(a +h) has changed from 4 to 2. 


Exercise 5.4.3 (Computational) Load the program INTEGR and eval- 
uate ik exp(—z”) dz using the options Gaussian 2 point and Gaussian 4 
point. Notice that at each stage you will be asked to specify the number 
of panels, where the total number of points used will be equal to the 
number of panels times the number of points in the formula. Compare 
the results obtained with those in exercise 5.4.2. 


Exercise 5.4.4 (Computational) Use the program INTEGR to estimate 


the integral 
iil 
— dz. 
I Jaz 


Why do the trapezium rule and Simpson’s rule fail? Compare the results 
produced by the open Newton—Cotes three-point rule and the Gauss-— 
Legendre methods. Which is the most efficient? 


> 5.5 Romberg integration 


Whenever we approximate an integral by a quadrature formula it is 
essential that some estimate of the quadrature error is made. We have 
seen that using the composite trapezium rule it is possible to show that 


b 
poms / f(x) dx = T(h) — 3h?(b- a) f"(8), aekn <b, ~ (5:51) 
where 
T(h) = $h[f(a) + 2f(a +h) +...+ f(d)]. 
If we halve h then we have 


I=T (5h) - (5h)? (6 - a) f"(0). (5.5.2) 
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Now multiply equation (5.5.2) by four and subtract (5.5.1) from the 
result. We obtain 


31 = AT(4h) — T(h) — h(a) F"(0") — F"(9)]. (5.5.3) 
If 0’ = @ then 
5 (47 (5h) — T(h)) = T (Zh) + 3 (T (Sh) —T(A)) 


should be a better approximation for J than either T(h) or T(#). This 
process is called Richardson’s method. 


Example 5.5.1 (Computational) 

Consider the integral ir exp(z?) dz. Using the program INTEGR we 
obtain values for the trapezium rule in table 5.5.1. Let us now apply 
Richardson’s method to the results in the second column. For example, 
applying equation (5.5.3) with h = 0.25 we have 


T(sh) + 3 (T(4h)-T(A)] 
1.469712276 — 1.490678861 


= 1.46971276 + a fern aero end = 1.462723415. 


Likewise, we can apply this technique to the remaining entries in the 
second column of table 5.5.1. Notice that we can achieve 6 correct 
decimal places using only 32 points instead of the 1649 predicted in 
example 5.3.1. 


Table 5.5.1 Richardson’s method applied to fe exp(x?) dz. 


N h Trapezium Error Richardson Error 
4 0.25 1.490678861 2.8 x 107 


8 0.125 1.4697 12276.057.1 xi105° 5) 1462723415007.) 10s 
16 0.0625 1.464420310 1.8x 107% 1.462656321 4.6 x 107° 
32 0.03125 1.463094102 4.4.x 10-4 1.462652033 2.8 x 10-7 
64__ 0.015625, 1.462762347 1.1 x.10-* | 1462651762) 1k58 1004 


Evercise 5.5.1 (Analytical) If S(h) denotes the result of applying Simp- 
son’s rule with step A to estimate an integral J = ihe f(x) dz then 

h4 

I — S(h) = — 


130° — a) OM Gs 0.< b. (5.5.4) 
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Show that 


165(4h) — S(h 1py , [SGA) = Sth 
a 2) sigan EI (5.5.5) 


gives a better approximation for I. 


Richardson’s method can be applied in a similar way to any other 
quadrature rule but it can be developed to produce a systematic ex- 
trapolation procedure. Let us assume that [(h) is a result produced by 
some quadrature rule when applied to an integral J and that the error 
associated with this approximation is given by 


I= I(h) + Ah? + O(h?*") (5.5.6) 
where A is some constant independent of h. Similarly, 
ies (sh) + A(zh)P + O(($A)?T"). (5.5.7) 


Now we multiply (5.5.7) by 2? and then subtract (5.5.6) from the result. 
This produces 


(2?—1)I = 2P?7(4h) — I(h) + O(AP*") (5.5.8) 
I = ge. (nett) (5.5.9) 


Therefore, if I(h) and I(5h) are of precision p— 1 then 


2P (Eh) — 1(h) 
2P — 1 


I(zh) ~ I(h) 


5.5.10 
RET ( ) 


— I(h) =f 
has degree of precision > p. This process can be repeated and is called 


Romberg extrapolation. 


Example §.5.2 (Analytical) 
It can be shown that if T(h) is the result of applying the trapezium rule 


to estimate an integral J = f f(x) dz then 
I =T(h)+ Ah? + Bh*+Ch* +... (5.5.11) 


where A, B, C, ..., are constants independent of h. These constants 
depend only upon f in the sense that A is related to f”, B to f*” and 
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Table 5.5.2 Romberg extrapolation applied to (e: exp(z”) dz. 


Nh T(h) __‘T(h, 3h)‘ T(h, $h, 4A) 
AR O20 1.490678861 
8 0.125 1.469712276 1.462723415 
16 0.0625 1.464420310 1.462656321 1.462651849 
32 0.03125 1.463094102 1.462652033 1.462651747 


so on. Therefore, provided f has a bounded second derivative we can 
eliminate the term Ah? to obtain equation (5.5.10) with p = 2, Le. 


T(h, 4h) = T(4h) + 3 [T($h) — T(A)], 
which has error proportional to h*. From this we have that 
T=T(h,th)+ BRA+ Ch’ +... 


where B’ is independent of h but does depend on the fourth derivative 
of f. If f has a bounded fourth derivative we can apply (5.5.10) to 
eliminate the term B’h* and so on. As an example let us use the program 
INTEGR to evaluate is exp(z”) dz. This particular integral is set up 
as default, therefore, select the option Romberg and start with h = 0.5 
to produce the Romberg tableau shown in table 5.5.2. Notice that we 
look down the columns for convergence because the terms across the 
rows must necessarily approach one another by the way in which they 
are defined. 


Although Romberg extrapolation can produce a dramatic increase in 
the accuracy of a quadrature rule it must be used with care. 


Example 5.5.3 (Analytical) 


Consider the integral if Vz log(x) dz. The integrand is shown in fig- 
ure 5.5.1 and is clearly defined at x = 0 but this function does not 
have a bounded second derivative throughout the interval and so the 
expression (5.5.11) is not valid. This integral can be carried out by 
making a change of variable, z = t?, which transforms the integral to 


fe 4t? log(t) dt, the integrand of which has bounded derivatives of all 
orders. 


ROMBERG INTEGRATION 137 
0.0 
QZ 
-0.4 
-0.6 
-0.8 
=) 


0.0 0.2 0.4 0.6 0.8 1.0 
Figure 5.5.1 The integrand \/z log(z) 


Romberg extrapolation is based on the existence of the expression 
(5.5.11), but even if such an expression is not valid it is sometimes 
possible to find an expression for the error. For example, for integrals 
of the form 


1 
ral x1/? lop(«)g(x) dz, 
0 


where g(z) is sufficiently smooth, the error in the trapezium rule, T(h), 
is 


I — T(h) = Ah? log(h) + Bh? + Ch? + Dh? log(h) + Eh? +... 


Given a sequence of values h;,h2,h3,... the terms in this series can be 
eliminated systematically. (See Fox and Mayers (1965).) 


Exercise 5.5.2 (Computational) Use the option Romberg in the program 
INTEGR to evaluate i z’ dz. Why do three extrapolations produce 
an exact answer? 


Erercise 5.5.3 (Analytical) If T(h) and S(h) are results produced by the 
trapezium and Simpson’s method show that 


S(4h) = T($h) + § [T($A) —T(A)] . 


Furthermore, show that the first extrapolation of Simpson’s rule, and 
therefore, the second extrapolation of the trapezium rule, produces the 
closed Newton—Cotes rule for four points given in table 5.2.1. (The 
resulting method is called Simpson’s Sth rule. ) 
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Exercise 5.5.4 (Computational) Use the program INTEGR with option 
Trapezium to estimate the integral fe exp(z?) dz. Repeat this exercise 
for the prompts Simpson's 1/3rd and Simpson's 3/8th. Confirm that 
these values are the first, second and third columns of the corresponding 
Romberg table. 


> 5.6 Adaptive integration 


With all quadrature formulae the choice of the step size, h, is critical. 
Too large a value and the result will be poor, too small a value and the 
amount of work required will be large and the result may be prone to 
rounding errors. In this section we will now outline an procedure which 
can automatically adjust the step size. To see the importance of this 
approach consider the following example. 


10 


O 2 4 6 8 10 


Figure 5.6.1 The integrand [(x — 5)? + 0.1]7?. 


Example 5.6.1 (Analytical) 
Let us estimate the integral 


10 1 
/ Goa 


The integrand is shown in figure 5.6.1. If we were to apply Simpson’s 
rule with a uniform step size h to produce an estimate S(h) then the 
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error involved is proportional to h4f*”. Throughout most of the interval 
f** will be very small; only in the neighbourhood of z = 5 will it be 
significant. However, to ensure that we achieve an overall precision € we 
need to select h such that 


Te fae 
180 


180¢€ hs 
h< ( , ) 
= ee are 


in which case we will be forced to use a small value for h even though 
the fourth derivative is very small except near x = 5. In order to over- 
come this difficulty we will now permit the sub-division of the range of 
integration and use non-uniformly spaced abscissa. Let 1?(h) be an ap- 


<e€ 


proximation for the integral J = ie f(x) dx which has degree of precision 


Dp, ee. 
Sy a es OY (5.6.1) 


If we require a final accuracy € then we should select h so that 


, 1/p 
ns aa! 


but in most cases |f+1)|,,.x is unknown. In order to overcome this 
we sub-divide the range of integration (a,b) into (a,c) and (c,b), where 
c= $(a+6), then 


Ie Ta(hh) = AGA) fOV(61) 
Ie —1g(th) = A(2A)P™ fe+D(6,), 


where a < 6; <c andc < 62 < 6. Therefore, there will exist a value 6” 
in (a, 6) such that 


+1 
I 1¢($h) — 129A) = ALO (0). 
If we subtract equation (5.6.1) from this then 


1>(h) — 12(2h) — 12(4h) & APH fA) (9)(1 — 27?) 


140 INTEGRATION 


provided 6 = 6’. This last equation enables us to estimate ft (0), or 
more precisely 


>(h) — (15($h) + 13(5h 
AhP+! f(e+1) (9) es (13( ) a (S ))) 
Hence, if 
[a(h) — 153) — Ie(G4)| < (2? — De 


then 
I — (15(4h) + 18(4R))| < 


as required. If not we reconsider each sub-interval independently and 
use a criterion of te in each sub-interval. 


Example 5.6.2 (Computational) 
If we use the program INTEGR with the option Simpson 1/3rd to 
evaluate the integral 


10 1 
———_———_— dz 
/ (x — 5)? +0.1 


then 512 intervals are required to ensure that two successive estimates 
agree to eight decimal places. Using the option Gaussian 4 point Rule 
128 panels of four points are needed to achieve the same accuracy and 
so there is little to be gained in using this method. However, using 
the adaptive procedure with a relative error tolerance of 107° only 161 
function evaluations are required. 


> 5.7 The Clenshaw—Curtis method 


We have seen in the previous sections how it is possible to estimate 
an integral by approximating the integrand by a suitable polynomial, 
integrate the approximation and then estimate the error. In order to 
obtain the best polynomial approximation, of a given order, it would 
seem obvious to use the minimax approximation. However to do this 
we need to expend more effort in finding the minimax approximation 
than we would use with a poorer approximation of higher order. As a 
compromise we could use an approximation which is almost minimax 
but which can be found a great deal easily, for example a Chebyshev 


approximation. This observation is the basis of the Clenshaw—Curtis 
method. 
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Let us assume that we have a Chebyshev approximation for a func- 


tion f(x) which is given by 
f(z) = 5 col + c17\(z) + c2To(x) mie eae 
then the indefinite integral of f(z) is given by 


J f@a= boo [ Tode ter [Ti(z)dz +... 


However, 
[toa = praeae4Co=T(2)+0o 
[ra = ears de? + 0. = Ina) +c1, 


where Co, C and C} are constants of integration, and 


1 Tn41 es | 
Pa jdreae ite ie > 2: 
/ oe 5 (a4 a4 | n>? 


This last result can be obtained as follows: 
[oe Che = [eosin cos~'x) dz 
= J ~ccs(n8) sin ao, (7=cos 0) 


= -} [isinin + 1)6 —sin(n — 1)6] dé 


eel fata — In-1 
sent aie 


(5.7.1) 


(5.7.2) 


(5.7.3) 


We now recall that T2, is an even function and that To,4; 1s an odd 


function and then use these results to show that 


Toes, eds and Zan lts =—a()) 


Substituting these expressions into equation (5.7.2) we determine 


1 
PC UE Er — 
[ fee a 2 |S + ben (735) +50 (Fos a 


Co C2 C4 


- WN So eer papa 
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and if we knew the Chebyshev coefficients of the function f(x) we could 
find the required integral. Notice that the above series will converge 
much faster than the Chebyshev series for the integrand because of the 
term r? — 1 in the denominator. Unfortunately, the determination of 
the Chebyshev coefficients is time consuming and in general will involve 
integrals of the form 


Ch = : fe cos(n@) f (cos 6) dé (5.7.4) 


which may be as complicated as the original integral. However, for the 
moment let us assume that we apply the composite trapezium rule with 
N points to the integral (5.7.4) then we find approximate coefficients c/, 
as 


2 
a rr (S90 +91+---+9n-1+59N), (5.7.5) 


where g; = cos(kna/N) f(cos(ka/N)). Notice that in this form we can 
take advantage of the fact that when the number of points is doubled 
from N to 2N we only need to evaluate g at a further N points. It is 
possible to show that 


Ch = Cn + Con—-n + Conn + C4n—n + CAN4n (5.7.6) 


and if the coefficients of the Chebyshev series approach zero rapidly then 
the approximate coefficients should produce a good approximation. 


Example 5.7.1 (Computational) 
Consider the integral 
/ exp(x”) dx 

=1 
for which the results in table 5.7.1 were obtained. From table 5.7.1 we 
see that provided N > n we obtain a very good approximation for the 
real coefficient which is given in the final column. If N = n then the 
approximation is out by a factor of 2, which can be confirmed by putting 
N =n in (5.7.6). We can now use these results to estimate the required 
integral. For example, if we include only the first four terms of the series 
given by equation (5.7.2) then 


1 
xp(22)dnw (a2 ee gee =5 UA 
[ exe?) ax EC aeeTs a 92532116, 


where the first neglected term is cg/63 % 1.746 x 107°. The exact value 
of this integral is 2.9253035 and the above estimate has error 1.76x 1075. 
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Table 5.7.1 The Chebyshev coefficients of exp(z”). 

N 2 4 8 16 Exact 
n=0 38.7182818 3.5078622 3.5067753 3.5067753 3.5067753 
n=2 1.7182818 0.8591409 0.8503917 0.8503917 0.8503917 
n=4 3.7182818 0.2104196 0.1052098 0.1052087 0.1052087 
n=6 1.7182818 0.8591410 0.0087492 0.0087221 0.0087221 
n=8 3.7182818 3.5078622 0.0010869 0.0005434 0.0005434 


Ezercise 5.7.1 (Computational) Use the option Trapezium in the pro- 
gram INTEGR to determine the first six Chebyshev coefficients for 
the function exp(—z”). How many function evaluations are required to 
determine the integral 
1 
/ exp(—z”) dz 

=i 
correct to four decimal places using the method of Clenshaw and Curtis? 
If the trapezium rule were applied directly to this integral how many in- 
tervals are required to determine a result which is correct to four decimal 
places? Which method is the most efficient and why? 


> Chapter 6 


> Differential Equations: Boundary 
Value Problems 


‘Pinning it down at both ends?’ 


> 6.1 Introduction 


In this chapter we shall discuss ways in which it is possible to approxi- 
mate the solution of differential equations of the form 


y = f(z,y,y'). (6.1.1) 


The general solution of this equation involves exactly two arbitrary con- 
stants. If two additional conditions are supplied then we may be able to 
find a unique solution, however, it is possible that the conditions which 
are supplied are not satisfied by any member of the class of functions de- 
termined by the general solution. If both conditions are specified at the 
same point then the problem is said to be an initial value problem and in 
this case under relatively simple conditions it is possible to tell whether 
or not a particular problem has a unique solution. If the two conditions 
are specified at different points then the problem is called a boundary 
value problem. There are relatively few cases when it is possible to say 
that a particular boundary value problem has a unique solution. (For a 
simple introduction see Harding and Quinney (1986).) Because of the 
difficulties involved with boundary value problems we will concentrate 
on linear problems which may be written in the form 


y" + p(z)y' + 4(2)y = r(z), (6.1.2) 
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ya)=a, y(b)=8, ards. 
Without loss of generality we may take a = 0 and b = 1 for if this is not 
the case we can always transform any finite interval (a,b) to (0,1) by 
means of a suitable mapping. We will also use the interval (—1,1) when 
it is more convenient to do so. 

In A Simple Introduction to Numerical Analysis it was shown how 
it was possible to determine an approximate solution of the two-point 
boundary value problem (6.1.2) by either a shooting method or by setting 
up a system of linear equations produced by replacing the derivatives in 
the differential equation by finite differences. In either case we obtain a 
representation of the solution only at a finite set of points; if we require 
the solution at a non-mesh point then we have to use some form of 
interpolation. In this chapter we will see how it is possible to obtain a 


solution in the form 
t= VN 


yn(z) = »S ci bi(Z), (6.1.3) 
2=0 
where each of the functions $;(zx) is defined for all points in the interval 
in which the solution is required. 


> 6.2 Solution in series 


We begin by considering equations of the form 
y'+p(z)y+a(z)y=r(z), O<2<1, (6.2.1) 


where p and q are rational functions. If p and q are continuous functions 
at a point ro then zo is said to be an ordinary point. If p and q are 
not continuous at zo, but (x — x9) p(x) and (x — 29)?q(z) are continuous 
at zo then zo is said to be a regular singularity. All other points are 
called irregular singularities. Near an ordinary point the solution of the 
differential equation can be expressed as a Taylor series expansion which 
converges as far as the next singularity. Near a regular singularity the 
solutions y; and y2 are of the form 


and 
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where 1, p2 and C, together with c, and dy, are to be determined. In 
both cases the series converges as far as the next discontinuity of either p 
or q. If C #0 then these solutions are linearly independent. However, it 
is possible to show that in all cases the general solution of the differential 
equation is given by 


y(x) = c1yi(2) + coy2(z). 


Since such solutions are known to exist we can substitute these expres- 
sions into the differential equation and try to determine values for p1, 
p2, C, and the coefficients c, and d,, n = 0,1,.... This is usually called 
Frobenius’s method or the method of solution in series. (See Bender and 
Orszag (1978).) 


Example 6.2.1 (Analytical) 
The differential equation 


ay" +2y=0 


has a regular singularity at ro = 0, therefore, we look for solutions of 
the form 


n=0o 
TGs) ee S Char 
n=0 
Substituting this expression into the differential equation we find that 
z?y" + 2y = co(p(p — 1) + 2)2? +.¢,(p? +p +2)2°t! +... 


However, the right-hand side must vanish identically, and so we select 
the coefficients and p so that this is the case. The coefficient of zr? is 


co(p’ — p+ 2) 


which is zero for all values of co provided p = 1 or p = 2. The coefficient 
of x°t? is c:(p? + p+ 1) and since p? + p41 # 0 for either p; or 
p = 2 we must have c; = 0. Furthermore, for these choices of p all the 
other coefficients are zero and so we obtain two solutions y;(z) = z and 
y2(z) = x”. The general solution is then given by 


y(z) = Az + Bz?, 


where A and B are arbitrary constants. Notice that even though the dif- 
ferential equation has a regular singularity at x = 0 the general solution 
is a continuous function near this point. 
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Erercise 6.2.1 (Analytical) Use the method of Frobenius to determine 
solutions of 


zy" —ay +y=0, 
which are valid for all values |z| < oo. 


Example 6.2.1 is a particularly simple example and in general it 
will be far more difficult to derive a simple expression for the solution. 
However, it does suggest that whenever we are looking for a continuous 
solution we should consider possible solutions in terms of a simple Taylor 
series expansion, 1.e. solutions of the form 


y(z) = 5 ene”. (6.2.2) 


n=0 


This type of expansion was considered in Chapter 3 and we have seen 
that in most cases it is more convenient to consider a Chebyshev expan- 
sion instead, therefore, we now look for solutions of the form 


nm=Co 


Gels SS Textnseah (6.2.3) 


n=0 


If we were to substitute equation (6.2.3) directly into the differential 


equation 
y + p(a)y’ +a(a)y.= r(x), (6.2.4) 


then we would need to differentiate each Chebyshev polynomial which 
requires considerable effort. Instead we integrate the differential equa- 
tion (6.2.4) to give 


ult py + [p+ quae = pravta 
and then again to produce 


yt | ovde [ [(-v + auaede = f [rdede+ Av +B, 


where A and B are arbitrary constants. We can now use some more 
properties of Chebyshev series to simplify this expression. In particular 


[oe dz = TEA — ——— n> 2 (6.2.5) 
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with 


and 


[ti@ae= fede = $2 = Pa) + Tole) 


(See equation (5.7.3).) Also 
1 =! A 
oT.) ( ; ) pee woe) (6.2.6) 


This last expression appears very complicated but it provides a key to 
determining the solution of the differential equation (6.2.4). 


Exercise 6.2.1 (Analytical) Deduce the expression (6.2.6). (Hint: recall 
that the Chebyshev polynomials are defined by 


Lie) S20 lee ener) 
which may be rearranged to give 
2Tn(z) = 3(Tr4(z) + Ta-1(2)) 


and this is equation (6.2.6) with r= 1. Now write z"T, = x’~1!(z7,).) 


Example 6.2.2 (Analytical/Computational) 
Consider the solution of the differential equation 


Yay te eee 


subject to the conditions y(0) = 1, y(1) = 3. Integrating the differential 
equation twice gives 


Substituting equation (6.2.3) the integral contains terms f{7T;, dz which 
can be replaced using equation (6.2.5). Next we replace the polyno- 
mial on the right-hand side of the differential equation by the equivalent 
Chebyshev series. For example, 


Teas Det Ta, 
2 pe ees loge 1): 


Dlr wle 
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Finally, we collect together all the coefficients of each Chebyshev poly- 
nomial, TJ, for n > 2. The coefficients of Tp and T; involve the arbitrary 
constants A and B and are therefore omitted. This gives 


Band Coyere ree) ee er 
cat (Sx6 — exe) = a 
at (SG - G8) = 0 
cs + (S53 — So) =. 0 


At x = 0 we use the boundary condition y(0) = 1 with 


0, n odd, 
EG) { (—1)"/?, n even, 


and so 
y(0) = 5co —coteg4—e6 +... = 1. 


Similarly, at z = 1 
“sl y= scot ci tco+e3+...=3. 


Therefore, we have an infinite set of linear equations in the unknowns 
Co, C1, ..., but as we expect them to tend to zero rapidly let us take the 
equations which come from the boundary conditions and the first four 
which are derived from the differential equation. All other coefficients 
are assumed to vanish. The result is a system of six linear equations in 
Co, C1, C2, C3, C4 and cs. Solving the resulting equations by Gaussian 
elimination with partial pivoting gives 


Co=2, cy = 2.045904223, 
co=0, cz = —0.04649799, 
cg =0, cs = 0.000593589. 


The solution of this problem is y(z) = 1+2z+ mas which has the 
Chebyshev expansion 


y(x) = 1+ 2.04590797T; — 0.046498073 + 0.000593675 — 0.0000035777. 
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Example 6.2.2 demonstrates how it is possible to use Chebyshev poly- 
nomials to produce approximate solutions of a limited class of linear 
differential equations. This method can be quite advantageous if more 
general boundary conditions are supplied. 


Exercise 6.2.2 (Analytical/Computational) Use Chebyshev polynomials 
to solve the differential equation 


y’ +aey +ey=l+2e4z2’, 
subject to the boundary conditions 
y(0)=1, (0) + 2y(1) — y(—1) = -1. 


The boundary conditions here are very complicated but the solution of 
the problem is very similar to that in example 6.2.1. (Hint: terms such 
as z”T,,(z) can be replaced using equation (6.2.6) and integrals such as 
f T,(2) dz can be written as $ f(In41 + In-1) dz.) 


> 6.3 Finite element methods 


Both Frobenius’s method and the Chebyshev expansion assume an ex- 
pansion of the solution as an infinite series. We will now consider a 
similar approach based on finite expansions. For the sake of clarity we 
will consider only linear differential equation which are written in the 
self-adjoint form 


(p(z)y')' + a(z)y = f(z), (6.3.1) 


subject to the homogeneous boundary conditions y(0) = y(1) = 0. It is 
relatively easy to show that any linear two-point boundary value problem 
subject to known boundary values can be reduced to this form. 


Example 6.3.1 (Analytical) 
As a simple example let us reduce the equation 


gz" +2'/4+¢5z=2—-1, 2(0)=0, 2(1)=1, 


to the self-adjoint form (6.3.1). The first step is to reduce the boundary 
conditions to a homogeneous form. To do this set y(x) = z(z) — x for 
then y(0) = z(0) -0 = 0 and y(1) = z(1) —1 = 0 as required. Similarly, 
the differential equation becomes 


zy! +y 4+ 229y=-l+2-27%. 
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If we rewrite equation (6.3.1) as py’ + p'y’ + qy = f and compare coef- 
ficients then we see that p’ = 1 and so p = z. Similarly, g(x) = z° and 
f(z) = -1+2-— 2%, and so the given equation becomes 


(zy) + 2°y = -14+ 25-27%, 
subject to y(0) = y(1) = 0. 
Exercise 6.3.1 (Analytical) Reduce the problem 
Chetealee Hee eau Z(G) Ome e( 1) nl 
to the self-adjoint form (6.3.1). 


The idea of the finite element method is to select a set of linearly 
independent basis functions ¢,;(r), i = 1,...,.N, each of which satisfies 
the homogeneous boundary conditions 


$i (0) = ¢;(1) = (i) 


and then to look for an approximate solution of the differential equation 


of the form 
1=N 


i) Se ern (6.3.2) 
i 

Notice that Yy(x) is defined at every point of the interval (0,1) and 
automatically satisfies the boundary conditions. The problem now is 
to decide which basis function to use and then how to determine the 
coefficients c;, 7 = 1,..,N, in such a way as to make Yy(z) the best 
possible approximation for the solution, i.e. we wish to minimize the 
error between y(x) and Yn(z). 


Example 6.3.2 
Two possible sets of basis functions are 


ror (2 eeesIt( 302.) bat Lyon apple 0 <1 


OO eee eet oe toils lV, OTS es 1, 


The problem of determining suitable values for the coefficients c; may 
be tackled in a variety of different ways. 
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> 6.3.1 The method of collocation 


Differentiating equation (6.3.2) and substituting the resulting expression 
into equation (6.3.1) we obtain 


ci s— 1) i=N 
P(t) D) eidi'(2) + p'(2) Dd) cidi(z) + a(2) D | evds(2) = F(2) 


t=) 


or oN 
Sle) 4 (e) + p'(@)4i(z) + a(2)4s(a)]es = F(2). 


If Yy (x) were the exact solution then this expression would be satisfied 
identically, but, since we have only N parameters at our disposal, 1.e. 
cj, 1 = 1,...,N, we cannot expect it to be satisfied at more than NV 
points. We can select any N points, z;,k = 1,...,N, but once they are 
specified the solution can be determined by solving the linear equations 


t=N 
Yd [P(ee) bi! (re) + P'(ze)d;(ee) + 4( xe) bi(ee ler = f(we), k= 1,...,N, 
t=1 

(6.3.3) 
by Gaussian elimination with partial pivoting, if necessary. This is called 
a collocation method and the points r,, k = 1,...,N, are called the 


collocation points. 


Example 6.3.3 (Computational) 

In order to demonstrate the collocation method let us take the basis 
function ¢;(z) = sin(i7z), 7 = 1,2,3, and then use the collocation points 
zz = 0.25, 0.5, 0.75, to compute an approximate solution for the bound- 
ary value problem 


: =e a y(0) = y(1) = 0. 


The trial solution is given by 
Y3(z) = c) sin wz + cosin27z + c3sin 372, 
where the coefficients c,, cp and c3 satisfy the linear equations 


Ac=b, 
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with the components of A given by 


: 2 
ane = (-(in) at =) SIN17TL, 
and those of b by 


ae (2r, = 4) 


ih SS ee 
pe sey Tere 2” 


| ema DPT 


Solving these equations gives 
c; = 0.161463, co =0.016691, cz = 0.0045005. 


This trial solution is compared with the exact solution of this problem, 
y(z) = x(1— 2)(1+2)7}, in figure 6.3.1. 


0.18 
0.16 
0.14 
0.12 
0.10 
0.08 
0.06 
0.04 


0.02 
0.00 


0.0 0.2 0.4 0.6 0.8 1.0 

Figure 6.3.1 Trial solutions Y3 and Y4 for example 6.3.3. 
Exercise 6.3.2 (Computational) Load the program FEM for which the 
default problem is the solution of the two-point boundary value problem 


2 (2x — 4) 
 — —<——$_—- yj = = Dale. 
DS (Glee aCe y(0) = yQ) 
Use the option Edit/inspect ODE to confirm that this is the case and 
then inspect the other options which are available at this point of the 
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program. The default problem can be examined using the prompt Op- 
tion: display problem and you can exit from this part of the program 
using the prompt Option: quit. After quitting the Edit/inspect ODE 
the prompt will become Method: Collocation. When this is the case 
press RETURN and then accept the option Edit/inspect basis. There 
are several sets of basis functions which are available but for the moment 
when the prompts shows Select basis: SIN(n*PI*x) press RETURN 
and then set the number of basis functions to three. Use the prompt 
Option: plot to examine the default basis which is supplied. When you 
quit from this option the program will display the prompt Points data: 
quit. Press the SPACE bar to see the various ways in which the points 
can be specified. For the moment let us use the same points as in ex- 
ample 6.3.3, i.e. 0.25, 0.5 and 0.75. These points can be input using 
the keyboard or using the option Points data: standard grid. Then 
use the prompt Points data: quit in order to compute the approximate 
solution using the prompt Option: Calculate solution. The progress 
towards finding the solution is indicated displaying the elements of the 
matrices A and b as they are found. The final stage is to plot the so- 
lution using the prompt Option: plot solution. The result should be as 
in figure 6.3.1. 


Exercise 6.3.3 (Computational) Use the program FEM to determine an 
approximate solution of the problem in exercise 6.3.2 using the basis 
{¢i(z) = sinimz}, « = 1,2,3,4. Compare the results with those in 
exercise 6.3.2. Try different sets of collocation points and determine the 
best possible choice. 


Exercise 6.3.4 (Computational) Load the program FEM and the option 
Edit/inspect basis to input the basis functions {¢;(z) = z‘(1 — z)}, 
t= 1,2,3,4. (You will also need to provide ¢; and ¢/’.) Use the method 
of collocation to determine approximate solutions of the default problem. 
Suggest alternative basis functions and compare the results which are 
obtained. 


In general the matrix associated with the collocation method, (6.3.3), 
will be full whereas with the finite difference method it is only a tridiag- 
onal matrix. It is possible to obtain basis functions which give a banded 
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matrix. As an example, consider the function 
0, zr<-2, 
+((2+2)°], —2<2<-l, 
Beet) eed (Lee) | ee <r <0, 
pee ees et) |e icat el, 


1(2 — 2)3, leer, 


=) 


) Sh. 


This function is called a cubic B-spline since it is a piecewise cubic which 
has continuous first and second derivatives. It is shown in figure 6.3.2. 


B(x) 


-2 =| O 1 2 


Figure 6.3.2 A cubic B-spline centred about z = 0. 


In order to construct a set of basis function ¢;(x), 7 = 1,...,N, we 
begin by partitioning the interval (0,1) by a set of N equally spaced 
points, 0 < 2; = ih < 1, (N+1)h= 1. For each point we define the 
basis function B;(x) = B((x — 2;)/h). Notice that outside the interval 
(x;~2, 2:42) the function B;(x) is identically zero. In order to apply the 
collocation method we need to evaluate the B-splines, and their first 
two derivatives, at the points z;, 7 = 1,...,N, but this turns out to 
be particularly easy as shown in table 6.3.1. However, as each B-spline 
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involves five points it is necessary to modify the basis functions near 
xz =0 and z =1. The basis function with modifications are 


B,(x) — £Bo(z), — ih. 
diz) = 4 Biz), 1<i<N-1, (6.3.5) 
By (2) 1 BN a) med aN 
to ensure that the boundary conditions are satisfied. Substituting these 
expressions into equation (6.3.3) we obtain the approximate solution 


i=N 


yn(t) = Ss c;9;(2), (6.3.6) 


Seal 


where 


Ac = b> 


withoe! =" [c1, (.-,em|o 9 DY = (fen; 2M )- Only themthrec 
main diagonals of the matrix A are non-zero and are given by 


ma = ~ spa) + ap(e1) + Palen), 
ao. = p5p(t2) ~ FP'(22) + a(22), 
(pare — S23) +425), a hee Ge 
Os ipa splay) + ap(2;) + 2D, 9 = ee Vig 
aM,M = — a 5p(em) - P'(em) + = a(em), 
aM-1.M = wP(2m-1) + = pi(zm—1) + (204-1). 


Exercise 6.3.5 (Computational) Load the program FEM to find an ap- 
proximate solution of the differential equation 


y (2x — 4) 


(+z) (+e) 


(y)_—2 


subject to the boundary conditions y(0) = 0, y(1) = 0, using cubic 
B-splines given by (6.3.4) and five equally spaced collocation points. 
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Table 6.3.1 Values of B;, B; and BY at 25,241. 


Tj-1 x; Ti+ 
B;(2) 1/4 1 1/4 
Bi (x) 3/4h 0 —3/4h 
B" (2) 3/2h? —3/h? 3/2h? 


To do this accept the options Edit/inspect basis? and Select basis: 
cubic B-splines. Set the number of splines to five and then accept the 
prompt Option: plot. Now plot the five splines which form the basis. 
This should produce a display like figure 6.3.3. Notice that the basis 
functions at either end are modified so that they include only 0, 1 and 
the collocation points. When you quit from this option the program will 
display the prompt Option: Calculate solution; press RETURN and 
the solution will then be determined. The computed solution can now 
be plotted using Option: plot solution. Compare the solution with that 
obtained in exercise 6.3.3. What is the minimum number of splines to 
get a reasonable approximation for the solution? Why are at least three 
splines required? 


) VW V3 V2 2/3 V6 1 


Figure 6.3.3 The five basis functions for exercise 6.3.5. 


The choice of the collocation points, r,, k = 1,..., N, is as impor- 
tant as the choice of basis functions. We have suggested that they be 
evenly spaced between 0 and 1 but this is by no means essential. Other 
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possibilities include the Gaussian quadrature points or the roots of a 
modified Chebyshev polynomial of degree NV. 


Erercise 6.3.6 (Computational) Load the program FEM with the de- 
fault problem and use the method of collocation with the basis functions 
{sin(nwz)}, n = 1,2,3,4 and collocation points as the roots of 74(zx). 
Compare the result with that obtained using equally spaced points. Why 
is the solution using the roots of T4 so poor? 


> 6.3.2 Galerkin methods 


We now look at an alternative way of minimizing the error between the 
trial solution Yy and the exact solution y(z). Let Qn be the set of 
functions which are generated by all possible linear combinations of a 
given set of basis functions, ¢;(z), 2 = 1,...,N. Therefore, a function 
Yn 1s a member of {2x if and only if 


i=N 
Y¥v(z) = )5 cdi(z). (6.3.7) 
Oss 
For each Yy we can determine the residual 


r = (pY{.)' +4Yn —f. (6.3.8) 


In general r will not be a member of Qa; however, there will be a 
projection of r onto Qn. Figure 6.3.4 illustrates this situation for N = 2. 


Figure 6.3.4 The projection of r onto Qy. 


The projection of r onto Qy will depend on the coefficients c¢; in 
equation (6.3.7) but it will be zero if they are selected so that r is 
orthogonal to each basis element. This is the basis of the Galerkin 
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method. We select the coefficients of equation (6.3.7) so that 


1 
| [(pY¥n)' +9¥n —flée(z)de=0, k=1,...,n. (6.3.9) 


Substituting for Yy from equation (6.3.7) and interchanging the sum 
and the integral produces a system of N linear equations given by 


Ac Db: 
where b; = He (2) or) dzy b= 1 nsand 


Aki = i [(p(x)¢;)' + q(x) di]ox dex. 


Integrating by parts reduces this expression to 


1 1 
Gy; = - | pla)dedide + | q(x) Gide. (6.3.10) 


Notice that the matrix A is symmetric, however, only if the basis func- 
tions are particularly simple is it possible to express its entries in a 
closed form. In many cases it will be necessary to resort to numerical 
integration to find the elements of A. Furthermore, in most cases the 
matrix A will be full, but there do exist functions for which it reduces 
to a tridiagonal form. 


Example 6.3.4 
Given a set of points x; = jh,j =0,...,N, Nh=1, the functions 
$(e@=2j-1), eja1 Se <az, 
)= 4 g(2j41-2), 2) SBS z41, (6.3.11) 
0, otherwise, 

are called linear B-splines and produce a tridiagonal system of linear 
equations. 
Exercise 6.3.7 (Analytical) Show that if the linear B-splines given in 
Example 6.3.5 are applied to the differential equation y” = f(x), subject 
to the boundary conditions y(0) = y(1) = 0, then the matrix A given in 
equation (6.3.10) reduces to 


—2 1 
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Exercise 6.3.3 shows the similarity between the Galerkin method and 
finite difference methods (see A Simple Introduction to Numerical Anal- 
ysis chapter 5). It is possible to show that the difference between the 
computed solution and the exact solution is proportional to h? provided 
the solution has sufficient bounded derivatives. If cubic splines are used 
instead then the Galerkin method is fourth order but the linear equations 
which have to be solved have five elements in most of the rows. 


Erercise 6.3.8 (Computational) Load the program FEM and the prob- 
lem 
(cy) +2°y=—-1+2—-2%, y(O)-= 1) =O, 


Using linear B-splines compute solutions for four basis functions for the 
option Galerkin. What happens as the number of basis functions in- 
creases’? Repeat this exercise with cubic B-splines. 


Exercise 6.3.9 (Computational) Load the program FEM with the data 
file FE which contains the problem 


y” — 167*y = —417” sin(57z) 


subject to homogeneous boundary conditions. Compute solutions for 
this problem using the options Method: Collocation and Method: 
Galerkin using the basis functions {sin(z), sin(37z), sin(57r)}. Com- 
pare the computed solution with that obtained using (i) cubic B-splines 
and (ii) linear B-splines. Next the basis functions {sin(nmx)}, n = 1,2,3 
and suitable collocation points. Now add sin(47z) to the basis and 
recompute the solution. Finally, add sin(5a7xr) and repeat this proce- 
dure. What is significant for n > 5? Why is it better to use the basis 
{sin(nwz)} with n = 1,3,5 rather than n = 1,...,5? 


> 6.3.3 Rayleigh—Ritz methods 
Consider the differential equation 
(p(z)y’)’ + a(z)y = f(z), (6.3.12) 


subject to y(0) = y(1) = 0. It is possible to show that the solution of this 
problem is identical with the function which minimizes the functional 


I(y) = ii [p(2)(x/)? — q(2)y? + 2fy] de. (6.3.13) 
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This is called a variational principle (see Elsgolc (1963), Jordan and 
Smith (1977)). In general the functional 


1 
(yj / F(a,y,y') dz 


has stationary points which satisfy 


d 
Pr! 
Exercise 6.3.10 (Analytical) Show that the stationary values of the func- 
tional (6.3.13) are equivalent to the solutions of the differential equa- 


tion (6.3.12). (Hint: look at I(y + dy) — I(y).) 


Therefore, in order to solve the two-point boundary value problem 
given by (6.3.12) we need to determine a function which produces the 
smallest value for the functional J. However, this is as difficult as the 
criginal problem and so we shall only attempt to look amongst a subset 
of all possible functions. For example, given any set of basis functions 
¢i(z), i= 1,..., N we could consider only functions of the form 


and select the coefficients c; to minimize the functional. To do this we 
substitute Yy into (6.3.12) and then 


HewoGo if [p(2)(V4,)? — 9(2)(¥w)? — 2f(2)¥w] de 


i i=N 2 i=N A 
= [ [ne (Soa) =G(z) Sete) 


1=1 


Differentiating I(c,;,...,cn) with respect to c,, k = 1,...,N, will pro- 
duce a system of linear equations in the unknown coefficients. These 
equations can then be solved using Gaussian elimination as before. It 
is possible to show that if the same basis functions are used for the 
Galerkin and Rayleigh—-Ritz methods then the solutions are identical. 
(See Strang and Fix (1973).) 


162 DIFFERENTIAL EQUATIONS 


Ezercise 6.3.11 (Analytical) Show that if the linear B-splines given in 
example 6.3.4 are used in a variational approach to solve the differential 
equation y” = f(x), subject to y(0) = y(1) = 0, then the system of 
linear equations which results is identical to that in exercise 6.3.7. 


> 6.4 Non-linear problems 


It is easy to extend the idea of finite element methods to non-linear 
equations, but the solution of the resulting non-linear equations for the 
coefficients is far from simple. For example, consider the differential 
equation 


eof (aay teen a0) oa yl eas (6.4.1) 


Now consider a trial function of the form Yy = Sen c; ¢i(z) which we 
substitute into equation (6.4.1) to find 


$f ater<1(oFaa Sur) 


s=1 $=1 


If Yy were the exact solution then this expression would be satisfied 
identically but we only have N parameters at our disposal. We can now 
use collocation to ensure that the above equation is satisfied at the points 
zy, k=1,...,N. This results in a system of N non-linear equations for 
the coefficients c;. For the Galerkin method we ensure that the residual 
is orthogonal to each of the basis functions, i.e. 


jae f(2,Y¥w, VD) orde = 0. 


Integrating the first term by parts and applying the boundary conditions 
to each basis function we find 


t= 1 1 i=N t=N 
-S> (| #10 de) c= / ip oe le Nae ‘) ox daz, 
i=l 0 0 i=1 i=1 
fork = 1,..., N. If the basis functions are splines then it is possible to 


reduce the problem but in general it is beyond the scope of this book. 
Further details can be found in Strang and Fix (1973). 


> Appendix A 


> Mathematical Background 


In this appendix we shall list some of the background mathematical 
material which has been used in the text. 


> A.1 Taylor’s series 


Suppose that the power series 
n=0O0 
Tose 


converges in some interval, -R < z—a< R, then the sum of the series 
has a value for each value of x and so defines a function. Therefore, we 
may write 


f(x) = ao +a;(2—a)+a(z—a)?+.... a-R<zr<atR. 
If we set z = a it is clear that ap = f(a). If we now differentiate with 
respect to z and set x = a, then a, = f’(a). Repeating this procedure 
gives 
fa) 
ni 


It is possible to show that provided the series converges in some interval 
then these steps are legitimate and so 


 LO%C@) 
Soe 


n! 


Aan = 


(ois (2 —a)”. (Avleb) 


n=0 
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This expression is called the Taylor series expansion of the function f 
about the point z =a. 


> A.2 Taylor’s series with remainder 


The infinite Taylor series expansion (A.1.1) is difficult to manipulate 
and, furthermore, the function f may not have an infinite number of 
derivatives. To overcome these difficulties we may use the following 
result. 


Theorem. Taylor’s series with remainder. 


Suppose that f, f’, f”,..., f%, and f+) are all continuous in some 
interval containing z = a and zx = 6; then there is a number @ between 
a and 6 such that 


f(b) = f(a) + fi(a)(b—a) +402(b- a)? + HM (b — a)3+ 
oy LN — a)" + Rn. 
where the remainder term FR, 1s given by 


CAO ESO) 


Hie (n+1)! 


Notice that R, depends on both a and b but once they are determined 
6 is fixed. 


> A.3 The mean value theorem 


Consider a Taylor’s series expansion with n = 1 and b = z then 


f(z) = f(a) +(x — a) f'(9) 
f(z) — f(a) = (2 - a) f'(9) 


which is called the mean value theorem. 


> A.4 Rolle’s theorem 


Let f be a continuous function in an interval (a,b) such that f(a) = 0 
and f(b) = 0; then there exists a value @ between a and b such that 


f"(9) = 0. 
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Proof. From the mean value theorem we have that 


f(b) — f(a) = (b— a) f'(4), 
where @ lies between a and b and since f(a) = 0 and f(b) = 0 we have 


that f’(@) = 0 as required. 


> A.5 Mean value theorem for integrals 


Let f and g be continuous functions over an interval (a,b) and assume 
that g(x) is positive; then there exists a value 6 between a and 6 such 
that 


i “HONORE AO / TOME 


Proof. As g is a continuous function between z = a and z = b there 
exist constants m and M such that 


m< g(r) < M. 


Accordingly, 
mf(z) < g(z)f(x) < Mf(z) 
and 


m / joe | Honey | HOE: 


Therefore, there exists a value M between m and M such that 


i HOWE r ee / Howe 


Furthermore, since g is continuous there is a value # such that g(@) = M 
which concludes the proof. 


> Appendix B 


> Finite Difference Approximations 


In order to construct finite difference approximations for derivatives we 
consider the Taylor series expansion of a function y(z) about z = a, for 
example, 


h2 
y(a + h) = y(a) + hy'(a) + Sy"). 
We now re-arrange this expression to give 
y(a + h) — y(a) 
h 
and so we can use the expression 
y(a +h) — y(a) 
h 


as an approximation for y’(a). Furthermore, the error in such an ap- 
proximation is proportional to h provided the function y has a bounded 
derivative. We say that this approximation is O(h). Alternatively, we 
can consider the Taylor series expansions 


= v(a) + 5u"(0) 


/ h? / h3 
y(a +h) = y(a) + hy’(a) + rye (a) + ary (2) Sear 
and 
/ he Wt h3 ws 
Veh) — aE hy a) ae (a) — 314 (a) ee 


Subtracting, the second expansion from the first we have 
h3 
y(a +h) — y(a— h) = 2hy'(a) +2=y"(a) +... 
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and so 
y(a +h) — y(a~h) 
2h 
may be used to approximate y/ (a) with error O(h?). 


We can extend the above approach to construct approximations for 
higher derivatives as follows: 


h? h3 ht. 
y(a +h) = y(a) + hy'(a) + ay (a) + Th SPORT bust a 
and 


/ h? M/ h? Mt h4 
y(a— h) = y(a) — hy'(a) + Fy (a) — A oF a 


3 ies. $e 


Adding these two expressions together we obtain 


how 
y(a +h) + y(a—h) = 2y(a) + h?y"(a) + DY er 


from which 
y(a + h) — 2y(a) + y(a — h) 
12 


is an approximation for y’(a) with error O(h?2). 
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This series of computer illustrated texts covers mathematics, 
science and engineering courses at advanced secondary school, 
college and university levels. 


A computer illustrated text is a textbook with integrated software. 
In a CIT the normal flow of explanation using conventional 
printed text is reinforced by computer software which can 
illustrate concepts both literally through graphics and figuratively 
by example. Computer graphics are very good at conveying a 
general qualitative understanding of a subject because they are 
much more flexible than static printed diagrams and allow 
greater reader involvement. 


The text is not specific to any particular microcomputer. Software 
is available for the BBC series of machines (40/80 track disc 
formats) and the IBM PC and compatible machines. 


The software is not available from your local bookshop. but is 
easily obtainable using the order form on page i. 


Approximation techniques are widely used in mathematics and 
applied physics, as exact solutions are frequently impossible to 
obtain. This companion to A Simple Introduction to Numerical 
Analysis extends the previous text to consider problems in 
interpolation and approximation. Topics covered include the 
construction of interpolating functions, the determination of 
polynomial: and rational function approximations, numerical 
quadrature and the solution of boundary value problems in 
ordinary differential equations. As with A Simple Introduction to 
Numerical Analysis the text is integrated with a software package 
which allows the reader to work through numerous examples. It 
is also possible to use the software to consider problems which 
are beyond the scope of the text. 


The authors’ expertise in combining text and software has 
resulted in a very readable work. 


A Simple Introduction to Numerical Analysis 2 is aimed at first- 
year students and advanced school students, but will be access- 
ible to anyone with a knowledge of mathematical methods. 
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